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Abstract. 


The  nature  of  interaction  in  multivariable  control  systems  is  examined  The 
concept  of  matrix  dominance  is  introduced  as  a  weaker  condition  than  diagonal 
dominance  for  the  application  of  the  multivariable  Nyquist  stability  theorems.  The 

Nyquist  exact  loci  design  technique,  which  uses  the  exact  transfer  function 
between  the  i'th  input-output  pair  when  the  i'th  loop  is  open  and  all  other 
loops  are  closed,  is  developed  and  applied  to  the  design  of  constant 

compensator/controllers  for  three  chemical  plant  models. 

Transmittances  in  a  closed  loop  multivariable  control  system  are  shown  to 
be  either  direct,  parallel,  interaction  or  disturbance  transmittances.  It  is  further 
shown,  that  interaction  transmittances  in  an  m-input,  m-output  multivariable 

system  can  be  expressed  as  a  sum  of  1st,  2nd,  up  to  m-1'th  order  terms  A 

critical  review  of  published  interaction  measures  shows,  that  the  relative  gain 

array  and  similar  indicators  do  not  measure  closed  loop  interaction.  The  direct 
Nyquist  array,  which  is  an  integral  part  of  a  multivariable  control  system  design 

technique,  is  shown  to  give  information  about  both  open  and  closed  loop 
interactions. 

The  property  of  matrix  dominance,  based  on  the  mathematical  theory  of 

M-matrices,  is  introduced  together  with  simple  numerical  and  graphical  tests  as 

a  weaker  condition  than  diagonal  dominance  for  application  of  the  multivariable 

Nyquist  stability  theorems.  Row  and  column  dominance  are  shown  to  be  dual 
with  respect  to  a  diagonal  similarity  transformation,  and  an  algorithm  is 
presented  for  transfering  dominance  between  rows  or  columns.  This  algorithm 
can  be  used  to  calculate  minimal  Gershgorin  radii. 

The  exact  transfer  function  h;(s)  between  the  i'th  input  and  the  i’th  output 
of  a  multivariable  system  when  the  i'th  loop  is  open  and  all  other  loops  are 

closed  is  derived  The  Nyquist  exact  loci  design  technique,  which  uses  polar 

plots  of  kj(s)h;(s),  is  developed  as  a  parallel  to  familiar  single  loop  design 
techniques.  It  is  shown,  that  stability  can  be  analysed  based  on  the  Nyquist 

exact  loci  provided  the  return  difference  matrix  is  matrix  dominant.  The  Nyquist 
exact  loci  give  direct  information  about  actual  bandwidth,  rise  time,  overshoot 


IV 


and  settling  time. 

The  Nyquist  exact  loci  (NEL)  design  technique  was  used  to  design 
constant  compensator/controllers  for  a  double  effect  evaporator,  an  open  loop 
unstable  chemical  reactor,  and  a  distillation  column  with  significant  time  delays. 
The  compensator/controllers  designed  using  the  NEL  procedure  are  equivalent  to 
those  designed  using  other  techniques. 
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1.  Introduction  and  objectives. 


1.1  Introduction. 

The  last  fifteen  years  have  brought  extensions  of  the  classical  single  va¬ 
riable  frequency  domain  control  system  design  techniques  of  Nyquist  and  Bode 
to  the  design  of  multivariable  control  systems.  These  developments  have  been 

spawned  by  practical  difficulties  in  the  use  of  modern  control  theory,  which  is 
to  a  large  extent  based  on  state  variable  formulations,  in  the  design  of  multiva¬ 
riable  process  control  systems,  Foss  (1973),  Horowitz  and  Shaked  (1975),  and 
Lee  and  Weekman  (1976). 

The  development  of  frequency  domain  methods  in  multivariable  control 
system  design  has  been  conveniently  summarized  in  a  publication  edited  by 

MacFarlane  (1979).  Two  well  known  multivariable  design  procedures,  the  inverse 
Nyquist  array  (INA)  method  due  to  Rosenbrock  (1969)  and  the  characteristic 
locus  (CL)  method  due  to  MacFarlane  and  Belletrutti  (1973)  have  been  applied  to 
industrial  problems,  but  their  use  has  been  somewhat  restricted  because  the 
requirement  of  diagonal  dominance  in  the  former  is  too  strict  for  many 
industrial  systems,  and  the  characteristic  values  or  eigenvalues  of  the  latter  have 
little  physical  meaning  to  the  industrial  designer.  Furthermore  the  non-diagonal 
controller/compensator  matrices  these  techniques  often  produce  may  reduce  the 

integrity  of  the  control  system,  especially  toward  actuator  failures.  Thus  there  is 
a  need  for  a  multivariable  control  system  design  method,  which  would  be  an 
exact  parallel  to  single  variable  frequency  domain  design,  and  in  a  meaningful 
way  use  familiar  terms  like  gain  margin,  phase  margin  and  crossover  frequency 
in  a  multivariable  context.  Also  the  method  should  allow  the  design  of  multiloop 
control  systems  for  most  chemical  processes,  thus  giving  high  integrity  toward 
both  actuator  and  transducer  failures.  In  the  Department  of  Chemical  Engineering 
at  the  University  of  Alberta  research  in  multivariable  frequency  domain  design 

techniques  was  initiated  by  Dr.  D  G.  Fisher  after  he  had  worked  with  Dr.  H.H. 
Rosenbrock  and  Dr.  AG.J.  MacFarlane  during  a  sabbatical  leave.  Kuon  (1975) 
conducted  a  survey  of  multivariable  frequency  domain  design  methods  and 
formally  developed  the  direct  Nyquist  array  technique  as  a  straightforward 
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design  method  for  multivariable  control  systems  that  to  the  practitioner  appears 
as  a  direct  extension  of  familiar  single  variable  techniques. 

In  parallel  with  the  development  of  multivariable  frequency  domain  design 
techniques  a  large  research  effort  has  gone  into  the  analysis  and  measurement 
of  interactions  in  multivariable  systems,  Rijnsdorp  (1965),  Bristol  (1966,  1968, 

1977,  1978,  1979)  and  others.  The  main  goal  of  this  research  has  been  the 

development  of  methods  for  selecting  the  best  pairing  of  manipulated  and 
controlled  variables  in  the  design  of  a  multiloop  control  system.  However,  these 
methods  provide  little  fundamental  insight  into  the  nature  of  interactions  and  the 
way  they  are  transmitted  through  a  multivariable  control  system.  Also,  even 

though  some  of  the  variable  pairing  methods  have  been  applied  successfully, 

they  do  give  incorrect  answers  for  certain  systems.  Therefore  the  need  exists 

for  a  better  practical  understanding  of  interactions  and  methods  for  reliably 

assessing  interaction  during  control  system  design. 


1.2  Objectives  of  this  research. 

Based  on  the  above  considerations  the  general  area  of  interest  in  this 

investigation  is  the  design  of  linear,  multivariable,  time-invariant,  control  systems 
in  the  frequency  domain.  The  main  transfer  function  elements  of  the  class  of 
systems  considered  are  shown  in  figure  1.1.  It  is  assumed  the  feedback 
transfer  function  matrix  (TFM)  H(s)  =  I,  the  identity  matrix,  since  this  represents 
no  loss  in  scope. 

For  the  class  of  systems  defined  by  figure  1.1  the  specific  goals  of 

this  research  are: 

i.  to  investigate  the  nature  of  closed  loop  signal  transmission  paths 

with  the  aim  of  clarifying  closed  loop  transmittances  and  interac¬ 

tions. 

ii.  to  investigate  the  use  of  the  direct  Nyquist  array  (DNA)  of  GP(s) 

as  a  tool  for  pairing  manipulated  and  controlled  variables,  and  the 

DNA  of  Q(s)  as  an  indicator  during  the  design  phase  of  closed 
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Figure 


ockdiagram  of  mam  transfer  function  elements >  in  £ i  rnumvanable 
>ntrol  system.  G„(s)  is  the  plant  TFM.  Gu(s!  is  the  '°^rT™ 
a  compensator  TFM,  and  Ka(s)  is  a  diagonal  controller  TFM.  Q(s) 
GP(s)*K1(s)  is  the  compensated  plant  TFM. 
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loop  interactions. 

iii.  to  review  published  measures  of  interactions,  especially  with 

respect  to  their  application  as  design  tools. 

iv.  to  formally  introduce  the  use  of  M-matrices,  as  a  TFM  property 

called  matrix  dominance,  in  conditions  for  multivariable  control 
system  stability  analysis  using  multivariable  Nyquist  array  techniques. 

v.  to  develop  graphical  displays  to  test  for  matrix  dominance  and  to 

aid  in  compensator  design  for  non-matrix  dominant  systems. 

vi.  to  develop  a  multivariable  control  system  design  technique,  which  is 

a  direct  extension  of  the  classical  single  variable  techniques  of 
Nyquist  and  Bode,  and  includes  meaningful  use  of  gain  margin, 

phase  margin,  crossover  frequency  etc.  in  the  multivariable  context. 

This  design  technique  will  use  the  transfer  function  hj(s),  which  is 
the  exact  multivariable  equivalent  of  the  single  variable  transfer 

function  g(s),  see  figure  1.2.  hj(s)  is  the  transfer  function  between 

the  i’th  input  and  the  i'th  output  when  loop  i  is  open  and  all  other 
loops  are  closed. 

vii.  to  apply  the  NEL  design  procedure  to  the  design  of  control 

systems  for  several  linear  plant  models,  and  evaluate  the  resultant 
controllers,  and  compare  the  results  with  those  of  other  design 

techniques. 

Although  it  is  not  included  as  part  of  this  thesis,  this  work  has  also  involved 
the  design  specification  for  an  interactive  computer  aided  control  system  design 
package,  Jensen  (1980).  Part  of  this  program  system  has  been  implemented,  and 
is  designed  to  assist  in  the  use  of  the  proposed  design  technique  and  other 
multivariable  frequency  domain  design  techniques,  such  as  the  inverse  Nyquist 
array  method,  the  characteristic  locus  method,  and  the  direct  Nyquist  array 
method.  The  programs  have  been  written  in  a  user  oriented  form,  and  the 
internal  structure  of  the  package  is  transparent  to  the  user. 
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Fiqure  1.2  Input-output  relationship  for  a  SISO  system  as  shown 
function  g(s),  and  one  pair  of  a  MIMO  system  as 
transfer  function  hj(s). 


►  y(s) 


►y^s) 


by  transfer 
shown  by 
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1.3  Structure  of  thesis. 

This  thesis  is  divided  into  six  chapters.  Chapter  two  contains  an  analysis 

of  interactions  in  closed  loop  linear  multivariable  control  systems.  The  relaxation 
of  the  diagonal  dominance  requirement  of  the  multivariable  Nyquist  array 
techniques  is  discussed  in  chapter  three.  Chapter  four  contains  the  development 
of  the  Nyquist  exact  locus  (NEL)  design  procedure.  In  chapter  five  the  proposed 
design  procedure  is  evaluated  by  application  to  three  different  systems. 

The  treatment  of  interaction,  dominance  and  the  NEL  method  in  chapters 

two,  three  and  four  are  written,  so  each  chapter  can  be  read  independently. 
This  means  there  is  some  overlap  or  redundancy  between  the  three  chapters. 

Conclusions  specific  to  each  topic  are  stated  at  the  end  of  each  chapter,  with 
the  overall  contributions  and  conclusions  given  in  chapter  six. 

The  application  of  the  NEL  method  in  chapter  five  is  written  to 
emphasize  the  practical  aspects  of  the  design  procedure  in  handling 
non-dominant,  open-loop  unstable  systems,  and  systems  with  pure  time  delays. 


'  *  1  .  a  ^  ; 


2.  Interaction  Analysis 


2.1  Introduction  and  scope. 

In  recent  years  a  great  deal  of  research  has  gone  into  the  analysis  and 

measurement  of  interactions  in  multivariable  systems,  Rijnsdorp  (1965),  Bristol 
(1966,  1968,  1977,  1978,  1979),  Davison  (1969,  1970),  Suchanti  and  Fournier 
(1973),  Witcher  and  McAvoy  (1977),  Tung  and  Edgar  (1977,  1978),  Kominek  and 
Smith  (1979),  Gagnepam  and  Seborg  (1979),  and  Jaaksoo  (1979).  The  main 
objective  of  this  research  was  to  develop  methods  for  choosing  the  best 
pairing  of  manipulated  and  controlled  variables  as  the  first  step  in  the  design  of 

a  multiloop  control  system.  However,  the  nature  of  interactions  has  not  been 
clearly  defined  nor  satisfactorily  incorporated  into  an  overall  control  system 
design  scheme.  Also  the  distinction  between  interactions  and  the  various 
input/output  transmission  paths  in  a  multivariable  system  has  not  been  clearly 
made. 

In  this  chapter  the  term  interaction  is  defined  and  the  nature  of  interac¬ 

tion  is  examined  with  respect  to  its  measurement,  implications  for  variable 
pairing,  and  its  use  in  a  multivariable  control  system  design  procedure.  A 
classification  of  transmittances  in  closed  loop  multivariable  control  systems  into 
direct,  parallel,  interaction  and  disturbance  transmittances  is  introduced.  Interaction 
transmittances  are  further  shown  to  be  a  combination  of  1st,  2nd  and  higher 
order  terms,  and  the  implications  of  low  interactions  for  disturbance  attenuation 

are  discussed.  A  link  between  interaction  transmittances  and  the  return  difference 
of  Bode  is  established 

The  direct  Nyquist  array  (DNA)  display  of  the  open  loop  transfer  function 
matrix  is  shown  to  give  a  measure  of  the  amount  of  interaction  to  be 

expected  in  the  closed  loop  system,  and  hence  to  be  a  useful  tool  for  pairing 
controlled  and  manipulated  variables. 

Finally  published  tools  for  measuring  interaction  are  critically  reviewed,  and 
it  is  shown,  that  several  of  them  do  not  measure  true  interaction.  Rijnsdorp's 

interaction  quotient  is  extended  to  be  applicable  to  systems  with  m-inputs  and 

m-outputs,  and  the  relation  between  Bristol's  relative  gain  array  and  the 
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multivariable  Nyquist  array  design  techniques  is  discussed,  and  demonstrated  by 
examples. 


2.2  The  nature  of  interaction. 

In  a  multivariable  system  one  output  is  in  general  influenced  by  more 

than  one  input  or  conversely  one  reference  input  (setpoint)  will  in  general 
influence  more  than  one  output.  This  is  normally  the  situation  both  in  the  open 
loop  and  in  the  closed  loop  system.  The  term  interaction  can  then  logically  be 

defined  as  follows  by  an  extension  of  MacFarlane  s  (1972)  definition: 

Def  i  nition:  Interactions  in  a  closed  loop  multivariable  system  are 
determined  by  the  transmittances  influencing  the  way  in  which  a 
reference  input  q(s)  affects  the  set  of  outputs  {y-(s):j£i},  or 

alternatively  the  transmittances  influencing  the  way  in  which  an 
output  y-(s)  is  affected  by  the  set  of  reference  inputs  £r-(s):j£i}. 

* J 

Although  interaction  arises  from  the  structure  of  the  open  loop  system,  it  is 
evident,  that  the  nature  of  interaction  is  best  analysed  by  considering  the  closed 

loop  transfer  function  matrix  of  a  multivariable  feedback  control  system,  such 

as  is  shown  in  figure  1.1.  As  noted  below  much  of  the  existing  literature  on 

interactions  consider  the  multivariable  system  with  the  i'th  loop  open.  This  does 
not  appear  to  have  much  physical  justification  in  interaction  analysis.  The  above 
definition  means,  that  any  transmittance  between  a  reference  input  rj(s)  and  the 
corresponding  output  yj(s)  is  not  an  interaction  transmittance. 

The  closed  loop  relationship  between  the  output  vector  y(s),  the 
reference  input  vector  r(s)  and  the  disturbance  vector  f  (s)  for  the  system  in 

figure  1.1  is 

y  =  (1+GpK,  K^H)”1  GPK1  Kar  +  (1+GpK*  K^Hf1  GLf 
=  (l+QK^H)**1  QK^  r  +  (l+QK^H)’1  GL| 

-  F^QK^r  +  F^Gj  =  Rr  +  RLf  (2.1) 

In  the  above  equation  and  those  following  the  Laplace  argument,  s,  is  omitted 

for  convenience.  From  equation  2.1  the  elements  r.-(s)  of  the  closed  loop 
transfer  function  matrix  R(s)  are  given  by 
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det  F 


(2.2) 


where  c(.(s)  is  the  (l,i)'th  cofactor  of  the  return  difference  matrix,  q^(s)  is  the 
(l,j)'th  element  of  Q(s),  and  kj(s)  is  the  j'th  diagonal  element  of  Ka(s).  A  similar 
expression  can  be  obtained  for  the  elements  of  the  matrix  Ru(s). 

The  i'th  row  of  the  closed  loop  transfer  function  matrix  R(s)  in  equation 
2. 1  can  be  viewed  as  a  single  loop  system  with  reference  input  rj(s),  other 

measurable  inputs  rj(s),j=  1,...,m,j£i,  disturbances  £(s),k=  1  ,...,p  and  output  yj(s). 
Furthermore  it  is  possible  to  express  the  diagonal  element  rj.(s)  of  R(s)  as  a 

sum  of  two  terms:  one  involving  only  kj(s)qfj(s)/det  F(s),  and  one  involving 

products  of  cofactors  of  the  i'th  column  of  F(s)  and  elements  of  the  i'th 

column  of  Q(s).  These  two  terms  can  be  viewed  as  two  signal  transmission 
paths  from  q(s)  to  yj(s),  one  being  influenced  only  by  the  diagonal  element 

kj(s)qjj(s)  of  Q(s)Ka(s),  and  one  being  influenced  by  two  or  more  offdiagonal 

elements  of  Q(s)Ka(s).  The  transmittances  along  the  two  paths  are  therefore 

labelled  respectively  the  direct  transmittance  and  the  parallel  transmittance.  Thus 
the  following  result  is  obtained: 

Remark  1:  Transmittances  in  a  closed  loop  multivariable  system  can 
be  classified  as  direct,  parallel,  interaction  and  disturbance  transmit¬ 
tance,  as  shown  schematically  in  figure  2.1  and  defined  mathemati¬ 
cally  below. 

The  i'th  output  can  be  expressed  as  a  sum  of  four  different  terms  due  to 

respectively  direct,  parallel,  interaction  and  disturbance  transmittances: 

k i,-  "L  =5ki  Mw 

y  =  - r,  +  E  -  r; 

det  F  k=  1  det  F 


+ 


m 


m 


j=1  k=  1 

j*i 


det  F  J  J 


p  m 

rns: 

i  is  —  i 


det  F^J 


(2.3) 


where  g  jK  is  the  (k,j)'th  element  of  GL(s).  Equation  2.3  shows,  that  parallel, 
interaction  and  disturbance  transmittances  are  all  dependent  on  the  properties  of 
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Figure  2.1  Blockdiagram  representation  to  show  classification  of  transmittances 
between  r*(s)  and  yj(s)  in  a  multivariable  system  (cf.  equation  2.3). 
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the  return  difference  matrix  F(s)  through  the  cofactors  c--(s)  of  this  matrix.  The 

J 

determinant  of  F(s),  det  F(s),  is  common  to  all  terms  and  can  thus  be  factored 
out.  Table  2.1  shows  the  elements  of  R(s)  for  a  2x2  and  a  3x3  system  with 
H(s)  =  I.  All  elements  are  given  with  the  determinant  of  the  return  difference 

matrix  as  common  denominator,  and  the  presence  of  direct,  parallel  and  interac¬ 
tion  transmittance  terms  are  evident.  From  equation  2. 1  it  can  be  seen,  that  as 

the  single  loop  gains  kj(s),  i=1,...,m  of  K^(s)  approach  infinity,  R(s)  approaches 

H(s)  ,  and  RL(s)  approaches  0.  This  leads  to  the  following  well  known  result: 
Remark  2:  Interaction  transmittances  and  disturbance  transmittances 
approach  zero  as  all  controller  gains  approach  infinity. 

The  significance  and  role  of  these  various  transmittance  terms  is  best 

appreciated  when  the  common  design  objectives  of  a  multivariable  control 
system  are  expressed  in  terms  of  them: 

Design  objectives:  The  compensator  (s)  and  the  feedback  con¬ 
trollers  {kj(s):i=1 . mj  are  normally  designed  to  meet  the  following 

objectives: 

i.  The  disturbance  transmittances  are  minimized.  (The  disturbance  trans¬ 

mittances  will  approach  zero  as  the  feedback  controller  gains 
approach  infinity). 

ii.  The  interaction  transmittances  are  minimized.  (The  interaction  trans¬ 
mittances  will  also  approach  zero  as  the  feedback  controller  gains 

approach  infinity). 

iii.  Ideally  the  parallel  transmittances  would  be  designed  to  complement 
the  direct  transmittances.  However,  they  cannot  be  designed 
independently  of  each  other  and  of  the  interaction  transmittances. 

iv.  The  direct  transmittances  are  designed  to  achieve  fast  servo 

control. 

Examination  of  table  2. 1  reveals,  that  the  numerators  of  the  interaction 

transmittances  r-(s),  i£j,  for  the  3x3  system  consist  of  two  terms,  which  are 

-  0 

respectively  1x1  and  2x2  minors  of  the  open  loop  transfer  function  matrix 

Q(s)K  (s).  The  1x1  minor  is  associated  with  interaction  transmitted  directly  from 

one  loop  to  another  loop  This  is  refered  to  as  a  first  order  interaction 
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Table  2.1  The  elements  of  the  closed  loop  transfer  function  matrix  R(s)  for  a 
2x2  and  a  3x3  system  with  open  loop  transfer  function  matrix 
Q(s)K^(s). 


2x2  system: 

rn  ~  (k1  q  1i  +  det  QK^)/det  F 
r,n  =  k»<Wdet  F 

rai  =  kid*i/det  F 
raa,  =  +  det  QK*l/det  F 

3x3  system: 

rn  =  <k-i«i«  t,.k!k?,q”  '  d»*cW  +  k.ks<cui  ^  + 

det  QKJ/det  F 


1Z 

=  ^A^\A 

+ 

kzk3*qiz  q33 

15 

~  *k3^13 

+ 

kzk3(qzzqi3 

Z1 

<r~ 

* 

cr 

T" 

II 

+ 

kik3(qziq33 

ZZ 

~  <kaqzz 
det  QK 

t/>ki(qazqii 
a)/det  F 

qi*qz3))/det  F 

f 

qi*qai>  +  k*k3<q**q33  -  q3*q*«>  + 


^^£3  +  k1  k3^  <1  ^-^3 

r3l  =  <M31  + 


ri2  =  <k3.qa2.  +  kik*<q«q 


*^33,  '  1  Z  H 11  H  33, 


r33  =  +  k3ki<q3.qn 

det  QKJ/det  F  ^ 


q^qi3))/det  F 

q^q«),/det  F 
q^i  qiz))/det  F 

q^q^  +  k3k^(q^^  ~  q^*1  + 


- 
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transmittance  The  2x2  minor  is  associated  with  the  transmittance  of  interaction 
from  one  loop  to  another  via  a  third  This  is  refered  to  as  a  second  order 
interaction  transmittance.  For  a  m'th  order  multivariable  system  the  interaction 
transmittances  will  contain  minors  of  the  open  loop  transfer  function  matrix  of 
order  up  to  m-1.  This  leads  to  the  following: 

Remark  3:  Interaction  transmittances  in  a  mxm  closed  loop  multiva¬ 
riable  system  can  be  expressed  as  a  sum  of  minors  of  1st,  2nd, 

and  (m-l)'th  order.  The  i'th  order  interaction  transmittances  are 
associated  with  ixi  minors  of  the  open  loop  transfer  function 
matrix. 

The  general  term  'associated  with’  is  used  deliberately  above,  since  it  is  not 

possible  to  associate  a  physical  transmission  path  with  just  the  numerators  of 
the  elements  of  R(s)  as  given  in  equation  2.2  or  table  2.1.  The  physical 

transmission  paths  are  revealed  by  dividing  the  numerators  by  det  F(s)  using 
long  division,  which  shows  there  exist  an  infinite  number  of  paths.  The  primary 

physical  transmission  paths  can  also  be  established  by  using  combinatorics  on  a 
signal  flow  graph.  Table  2.2  gives  the  transmittances  along  the  primary 

transmission  paths  for  a  2x2  and  a  3x3  stystem.  The  term  'primary'  is  meant 
to  indicate,  that  no  section  of  the  path  is  traversed  more  than  once.  The 

transmittances  given  in  table  2.1  also  show  terms,  which  could  be  classified  1st 
and  2nd  order  interactions.  Compared  to  the  classification  introduced  above  this 
would  give  an  infinite  number  of  interaction  terms,  which  would  make  analysis 
of  interaction  more  complicated. 

The  elements  of  R(s)  can  be  calculated  from  the  the  minors  of  the  open 

loop  transfer  function  matrix  Q(s)K^(s)  using  the  following  general  expression: 
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r..  = - 

det  F 


m 


oj  +  z  q]I; 

4  i,=i  J’ 


+ 


m 

'i  =  1  U 


m  ..  . 

T.  Q'VfV 

=1,  +  1^ 


m 

+z: 

I  4  —  1 


Q 


"1- 

4*1 

+ 1 


+  <$••  detQ' 

0 


(2.4) 


where  1^!^,...  are  always  different  from  i  and  j.  In  equation  2.4  Qj  denotes  the 
1x1  minor  of  Q(s)Ka(s)  formed  by  deleting  all  rows  except  row  i  and  all 


14 


Table  2.2  The  transmittances  along  the  primary  physical  transmission  paths  in  a 
2x2  and  a  3x3  closed  loop  system  with  open  loop  transfer 
function  matrix  Q(s)K^(s). 


2x2  system: 
r;  = 

riz  “  k«,qi2, 
rZ1  k  1^2,1 

rzz  ~  ^z^z-z  ~  ki  k^/c!ai  3  iz 
3x3  system: 


r11 

=  k,k3q31q,s 

k1  kACli1 

=  kaqla,  '  kj-ka^a^ia 

r4 

k3q,j  "  kAk 3q 

r*i 

=  k^Q^  “  k*k3q31q33 

=  kj.q*,!,  _  kAk,q1s_q31 
kikak3(q13qiaqai  + 

=  k3qai  -  k,  k3qai  q,j 

r31 

k^q31  "  k1k^>qaiq3a 

rzz 

=  kaq„  -  k,kaq31qu 

rU 

=  k3q3a  “  k^k^q^^q^^ 
ki  k3k3(q1;tqA3q31  + 

k3k3.q  3^32, 
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columns  except  column  j,  £;•  is  the  Kronecker  delta,  and  =1  for  i=j, 

M  vj 

=0  for  iAj. 

Much  of  the  interaction  analysis  literature  is  concerned  with  the  transmit- 
tances  in  a  multivariable  system  in  which  the  i'th  loop  is  open  and  all  other 
loops  are  closed.  In  this  situation  the  closed  loop  relationship  between  the 
reference  input  vector  r(s),  the  disturbance  vector  |  (s)  and  the  output  vector 
y(s)  is 

y  =  (l+QK^HSr1  QK^r  4  (PQK^HSf1  G.  I 

la  ^ 

=  RV  4  R'|  (2.5) 

where  S  is  a  diagonal  matrix  with  the  i'th  diagonal  element  zero  and  all  other 
diagonal  elements  unity.  If  k  (s)  =  1  in  equation  2.5,  then  the  element  rji(s)  of 
R'(s)  is  the  transfer  function  for  which  the  i'th  loop  controller  k*(s)  should  be 
designed  to  take  into  account  the  multivariable  nature  of  the  system.  A  design 
procedure,  which  uses  these  transfer  functions,  is  developed  in  chapter  four 
and  applied  to  several  systems  in  chapter  five  An  expression  similar  to 
equation  2.3  for  the  i'th  output  when  all  loops  except  the  i'th  are  closed  can 
be  stated  as  follows: 


Yj  =  9ukir; 


4 


Likir; 


4 


m 

j=1 

j*i 


r.-  r- 

•J  J 


4 


P  , 

pi  ^ 


(2.6) 


where  Lj(s)  is  the  parallel  transmittance  between  uj(s)  and  y.(s), 
k-(s)rj(s),  when  the  i'th  loop  is  open.  The  parallel  transmittance 
calculated  by  the  following  expression,  which  is  derived  in  appendix 

m  m 


Z  ^  q5yk.  =  Z 

j#i  j*i 


and 

L;(S) 

A: 


Uj(s)  = 
can  be 


(2.7) 


where  (s)  is  the  ratio  of  the  (i,j)'th  cofactor  to  the  (i,i)'th  cofactor  of  the 

return  difference  matrix  F(s).  Rosenbrock  (1972)  gives  an  expression  for  hm(s)  = 
q^^s)  +  L^Js),  which  is  the  transfer  function  with  the  m'th  loop  open  and  all 

other  loops  closed  From  this  an  expression  for  L^s)  can  easily  be  derived 
However,  equation  2.7  is  simpler  and  more  flexible,  since  any  one  loop  can  be 

open,  whereas  Rosenbrock's  result  requires  the  m'th  loop  to  be  open. 
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As  discussed  later  the  transmittances  h.(s),  i=1,...,m  in  the  limit  as  the 
gains  in  all  closed  loops  approach  infinity  plays  an  important  part  in  several 

published  measures  of  interaction.  However,  from  the  foregoing  treatment  it  is 

evident  these  limiting  transmittances  do  not  represent  interaction  transmittances, 
but,  as  defined  in  figure  2.1,  only  the  limiting  values  of  the  sum  of  direct  and 

parallel  transmittance  with  one  loop  open  and  the  rest  closed.  In  the  limit  as 
the  gain  of  all  closed  loops  approach  infinity  the  ratio  of  cofactors  ($>••  (s) 

J 

approach  a  corresponding  ratio  of  cofactors  of  the  compensated  plant  transfer 
function  matrix  Q(s).  This  means  the  parallel  transmittance  L.(s)  does  not  in 
general  approach  zero  under  high  gain  feedback  in  all  closed  loops,  and  hence 
should  not  be  used  as  a  measure  of  true  interaction 

Equations  2.3  and  2.6  also  show,  that  parallel  transmittance  may  prove 
helpful  by  increasing  the  controller  bandwidth,  whereas  interaction  transmittances 
must  always  be  considered  harmful.  Since  r.(s)  and  r.(s)  ,  j^i,  can  have  opposing 

effects  on  y.(s),  and  the  disturbances  are  in  general  arbitrary  and  unknown  it  is 

not  possible  to  establish  a  priori  whether  a  particular  interaction  transmittance  is 

helpful  or  harmful.  Ideally  one  would  like  to  take  advantage  of  the  parallel 
transmittance  and  at  the  same  time  minimize  interaction  transmittances  in  all 
control  loops.  However,  these  are  conflicting  objectives,  as  evident  by  inspection 
of  table  2.2.  In  practice  it  will  only  be  possible  to  take  advantage  of  the 

parallel  transmittance  in  the  design,  or  to  minimize  interaction  transmittances  in 
some  loops  at  the  expense  of  more  interaction  in  other  loops.  Tools  for 
accomplishing  these  objectives  are  developed  in  chapters  three  and  four. 

Although  the  above  treatment  has  neglected  disturbance  transmittances 
they  are  closely  related  to  interaction  transmittances  In  the  servo  problem  the 
objective  is  to  make  a  change  in  one  input  variable  affect  only  one  output  va¬ 
riable,  i.e.  minimize  all  interaction  transmittances.  In  the  regulator  problem  the 

objective  is  to  minimize  the  influence  of  a  disturbance  on  all  output  variables, 
i.e.  minimize  all  disturbance  transmittances.  Equation  2.3  shows,  that  the  cofactors 

c--(s)  of  the  return  difference  matrix  influence  the  interaction  and  disturbance 

’J 

transmittances  in  identical  ways.  Hence  the  objective  of  both  the  servo  and  the 

regulator  problem  can  be  met  by  minimizing  the  cofactors  of  the  return 
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difference  matrix.  This  would  also  minimize  the  parallel  transmittances,  but  a 
practical  procedure  for  affecting  the  minimization  has  not  been  developed  yet. 


2.3  The  direct  Nyquist  array  (DNA) 

In  the  multivariable  Nyquist  type  design  techniques,  such  as  the  direct 
Nyquist  array  (DNA)  technique  and  the  inverse  Nyquist  array  (INA)  technique, 

classical  single  loop  control  system  design  approaches  have  been  extended  to 
the  design  of  multivariable  control  systems.  The  use  of  the  DNA  display,  which 
is  an  integral  part  of  an  established  multivariable  design  technique,  to  give 
information  about  closed  loop  interaction  and  pairing  of  variables  is  discussed 
below. 

2.3.1  The  DNA  display  and  interaction. 

The  DNA  display  in  the  direct  Nyquist  array  technique  is  used  to  design 
a  compensator  K^(s),  which  reduces  the  interaction  in  the  open  loop  system. 
The  implication  is  that  reducing  open  loop  interaction  will  reduce  closed  loop 
interaction  as  well.  In  the  previous  section  it  was  shown,  that  interaction  trans¬ 
mittances  can  be  classified  as  1st,  2nd.  and  higher  order,  and  that  i'th  order 
interaction  transmittances  are  associated  with  ixi  minors  of  the  open  loop 

transfer  function  matrix.  For  most  practical  systems,  such  as  chemical  process 
systems,  the  magnitude  of  the  minors  will  normally  decrease  as  the  order 
increases.  Hence, 

Remark  4:  In  most  practical  systems  the  first  order  interaction 

transmittances  provide  a  reasonable  approximation  to  the  total  in¬ 
teraction,  i.e  the  higher  order  interaction  transmittances  are 
generally  negligible 

Since  only  first  order  interactions  are  found  in  2x2  systems  the  above  remark 
is  clearly  valid  for  all  2x2  systems.  It  is  however  possible  to  construct 

examples  of  3x3  systems  which  contradict  the  above  statement,  i.e.  have  IQ'] I  > 
IQ'jJl  for  some  (i,j,l),  l£i,  l/j.  It  is  nevertheless  always  possible  to  compare  direct 
transmittance  and  first  order  interaction  transmittances  in  the  closed  loop  system 
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by  comparing  a  diagonal  element  and  the  corresponding  offdiagonal  elements  of 
the  open  loop  transfer  function  matrix  Q(s)K^(s).  The  direct  Nyquist  array 

displays  the  elements  of  Q(s)K^(s)  as  a  function  of  frequency  in  the  form  of 
an  array  of  polar  plots,  hence 

Remark  5:  The  DNA  display  of  Q(s)  gives  an  exact  measure  of 
direct  and  interaction  transmittances  in  the  open  loop  system,  as 
well  as  an  indication  of  the  direct  transmittance  and  the  first  order 

interaction  transmittances  in  the  dosed  loop  system. 

Thus  the  DNA  display  of  Q(s)  gives  an  indication  of  the  closed  loop  interac¬ 
tions  The  DNA  plot  of  Q(s)  further  has  the  advantage  of  being  an  integral  part 

of  an  established  multivariable  control  system  design  technique.  No  additional 
numerical  calculations  are  necessary  to  get  the  first  order  interaction  information, 
and  the  change  in  the  amount  of  interaction  can  be  followed  as  the  control 

system  design  progresses.  If  desired  higher  order  interactions  can  be  calculated 
from  the  elements  of  Q(s)K^(s)  and  displayed  as  part  of  the  DNA  design 
procedure. 

The  INA  does  not  offer  as  simple  a  physical  interpretation  of  the  display 
due  to  the  complicated  relationship  between  the  elements  of  R(s)  and  of  R_1(s). 
Furthermore  Clement  (1980)  has  demonstrated,  that  the  DNA-design  procedure  is 
preferable,  when  the  plant  transfer  function  matrix  is  experimentally  determined. 

2.3.2  The  DNA  display  and  variable  pairing. 

The  first  step  in  the  frequency  domain  design  of  a  multivariable  control 
system  usually  is  the  pairing  of  manipulated  and  controlled  variables.  The  aim  is 
to  find  the  pairing,  which  makes  it  possible  to  realize  the  control  system 
objectives,  eg  minimize  interactions  and  maximize  disturbance  attenuation,  with 
the  minimum  amount  of  effort.  Since  the  DNA  display  gives  a  measure  of 
interaction  the  pairing  objectives  can  be  approached  as  follows: 

Remark  6:  Pair  manipulated  and  controlled  variables  so  the  largest 
elements  of  the  DNA  display  of  Gp(s)  lie  on  the  main  diagonal. 

This  method  will  maximize  the  direct  transmittance  and  minimize  the  first  order 
interaction  transmittances  of  the  closed  loop  system.  Also  the  above  pairing 


19 


procedure  is  consistent  with  Bristol's  (1977)  relative  dynamic  gain  array 
approach,  but  requires  less  numerical  calculation.  This  is  the  case,  because  Gp(s) 
being  matrix  dominant  (cf.  chapter  3)  implies,  that  the  relative  dynamic  gain  array 
is  matrix  dominant,  Fiedler  and  Ptak  (1966). 


2.4  Critical  review  of  published  interaction  measures. 

Most  published  measures  of  interaction  have  shortcomings  such  as:  i)  they 
are  based  on  a  system  with  one  loop  open  rather  than  a  fully  closed  loop 
system,  and  ii)  they  measure  parallel  transmittance  rather  than  the  true  closed 
loop  interaction. 


2.4.1  The  interaction  quotient. 

One  of  the  first  suggested  measures  of  process  interaction,  was  the  in¬ 
teraction  quotient  proposed  by  Rijnsdorp  (1965).  Rijnsdorp  considered  only  2x2 
systems  with  emphasis  on  distillation  column  control.  For  2x2  systems  the  inter¬ 
action  quotient  is  defined  to  be 


x  =  SaiS^/g,,  g«,  (2.8) 

where  g;j(s)  is  an  element  of  the  plant  transfer  function  matrix  Gp(s).  The 
quotient  is  evaluated  for  s  =  i  co  and  plotted  as  a  polar  plot.  The  interaction 
quotient  has  the  properties,  that  it  is  dimensionless  and  invariant  under  scaling. 
Rijnsdorp  concludes,  that  when  the  static  value  AC  ( 0)  is  close  to  unity  interac¬ 
tion  causes  poor  control.  This  is  also  the  case  when  /C(s)  shows  increasing 

negative  phase  with  frequency,  and  a  magnitude  close  to  one,  e.g  a  pure  time 

delay.  However,  when  AC  (s)  is  approximately  constant  and  has  a  negative  real 
part  good  control  can  be  achieved  with  a  multiloop  control  system.  Poor 
integrity  results  if  AC  (s)  is  approximately  constant  and  has  a  real  part  greater 

than  one. 

The  relationship,  assuming  K^(s)  =  I, 


K  - - -  lim- 

9 ii  3**. 


'Z1 


lim- 
k> 


”1  +ki  9  ii 


L. 


■=  lim- 


ki 


+ 


9n 


k^°°q^ 


(2.9) 


' 


20 


shows  the  interaction  quotient  is  the  limiting  value  of  the  ratio  of  parallel 

transmittance  to  direct  transmittance  for  a  system  with  one  loop  open  and  one 
loop  closed  Hence  the  interaction  quotient  is  not  a  measure  of  interaction,  but 
a  measure  of  the  significance  of  parallel  transmittance  with  one  loop  open. 

Rijnsdorp's  conclusions  about  the  interaction  quotient  provide  an  answer  to  the 
question  Are  problems  likely  to  be  encountered  in  the  design  of  a  multiloop 
control  system  for  this  plant?'.  However,  the  interaction  quotient  does  not  give 
a  reliable  answer  to  the  question  'Is  interaction  (as  defined  in  section  2.1)  signi¬ 
ficant?’.  This  will  be  further  demonstrated  by  examples  in  section  2.5. 

The  interaction  quotient  as  an  indicator  of  possible  difficulties  in  multiloop 

design  can  be  extended  to  mxm  system  by  the  following  definition: 

m  m 


pi  ^\iqM 


h  ^  q^ 

j* 


(2.10) 
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where  7|p(s)  is  the  ratio  of  the  (i,j)'th  cofactor  to  the  (i,i)'th  cofactor  of  the 

compensated  plant  transfer  function  matrix  Q(s).  The  above  extended  definition 
follows  from  the  expression  for  Lj(s)  and  the  fact,  that  as  k  in  Kx{s)  =  kl 

approaches  infinity  the  ratio  of  cofactors  of  the  return  difference  matrix 

approach  the  corresponding  ratio  of  cofactors  of  the  compensated  plant 

transfer  function  matrix,  i.e  with  K^(s)  =  kl 

lim  4V,  =  ’H’fi  (2.11) 

k-*~=>  J  0 


The  proof  of  the  above  relationship  is  given  in  appendix  B.  It  is  easily  verified, 

that  the  extended  interaction  quotient  defined  by  equation  2.10  reduces  to 
Rijnsdorp's  for  a  2x2  system.  The  extended  interaction  quotient  has  the  same 

properties  as  Rijnsdorp’s  measure,  and  in  addition  reordering  of  inputs  or 

outputs  corresponds  to  an  equivalent  reordering  of  the  set  {l-C-  :i=1,...m}.  The 
extended  interaction  quotient  also  is  only  an  indicator  of  possible  difficulties  in 

multiloop  control  system  design,  and  not  a  measure  of  closed  loop  interaction. 


* 
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It  gives  a  measure  for  each  row  or  column  of  the  compensated  plant  transfer 
function  matrix,  that  is,  it  gives  a  measure  for  each  control  loop. 

Kommek  and  Smith  (1978)  gives  the  first  extensive  dynamic  interpretation 
of  the  interaction  quotient  They  use  polar  plots  of  K  (s),  s  =  i  to  ,  with  a 
superimposed  unit  circle.  They  classify  the  plots  as  showing  favourable  and 
unfavourable  interactions.  From  the  preceding  treatment  it  is  evident,  that 
Kominek  and  Smith  are  actually  judging  whether  or  not  parallel  transmittance  will 
be  advantageous  or  disadvantageous  in  a  multiloop  control  system. 


2.4.2  The  relative  gain  array. 

Bristol  (1966,1967)  suggested  a  steady  state  interaction  measure,  which  is 
a  relative  gain  ratio: 

r  (  ^  y ;/  U : ) I  Aub  =  0  for  k*jl 

JU-  ;:=  r - 1 - 5 - 4  (2  1  21 

^  L  (  ^  yj/  4^  uj)l  ^  y^  =  0  for  k£i  J 


The  numerator  corresponds  to  the  open  loop  gain  between  u-  and  y-  and  the 
denominator  is  the  gain  between  uj  and  y  j  with  all  other  outputs  perfectly 
controlled.  M  =  W  is  called  the  relative  gain  array  (RGA).  The  RGA,  as 

Rijnsdorp's  interaction  quotient,  is  restricted  to  plants  or  compensated  plants 
with  an  equal  number  of  inputs  and  outputs.  By  virtue  of  the  defining  equation 
the  elements  of  the  relative  gain  array  are  dimensionless.  This  imply  they  are 
invariant  under  scaling,  and  since  selection  of  final  proportional  controller  gams 
is  a  scaling  operation,  the  relative  gain  array  is  independent  of  the  final 
multiloop  controller  gains  K^(s).  The  RGA  further  has  the  property,  that  any  row 
or  column  sums  to  one,  and  that  reordering  inputs  or  outputs  implies  a  similar 
reordering  of  columns  or  rows  of  M.  Bristol  has  shown,  that  only  the  open 
loop  gains  are  necessary  to  evaluate  the  RGA,  since 

M  =  Q(0)  o[Q(0f1]T  =  tqi:(0^q.j(0)}  (2.13) 


where  the  symbol  ®  indicates  the  Hadamard  or  Schur  product  of  the  two 
matrices,  see  Johnson  (1974)  Bristol  (1968)  recommends,  that  controlled  and 
manipulated  variables  with  /u.^  positive  and  close  to  unity  be  paired.  Also  if 
any  u..  is  much  larger  than  one  or  less  than  zero  pairing  the  corresponding 
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variables  will  result  in  a  loop,  which  is  difficult  to  control.  Bristol  (1977)  has 
also  shown,  that  pairings,  which  gives  negative  relative  gains  on  the  diagonal 
should  be  avoided,  since  this  would  result  in  a  control  loop  with  nonminimum 
phase  characteristics. 

The  steady  state  relative  gain  array  has  been  successfully  applied  by 
several  authors,  e.g.  Stainthorp  (1972),  Nisenfeld  (1973),  as  a  variable  pairing 
tool  in  designing  multiloop  systems  to  give  minimum  closed  loop  interaction. 

The  steady  state  nature  of  the  RGA  neglects  any  high  frequency 
dynamics,  which  could  be  important  in  certain  systems,  such  as  the  turbo 
alternator  considered  by  Ahson  and  Nicholson  (1976)  and  by  Taiwo  (1978)  and 
used  in  example  3  of  section  2.5.  Witcher  and  McAvoy  (1977)  and  later  Bristol 
(1978)  therefore  suggested  an  intuitive  dynamic  extension  by  defining  a  relative 
dynamic  gain  array  (RDGA)  as  follows 

M(s)  =  Q(s)  o  [Q(s)-1I1T  =  {qj-(s)*cjjj(s)}  (2.14) 

which  posseses  the  same  properties  as  the  steady  state  equivalent,  e.g. 
dimensionless.  Witcher  and  McAvoy  (1977)  state,  that  interaction  is  small  if  the 
magnitude  of  the  diagonal  elements  of  M(s)  are  close  to  one  and  the  magnitude 
of  all  other  elements  are  small.  This  is  equivalent  to  saying  M(s)  should  be  a 
diagonally  dominant  matrix. 

The  relationship,  for  a  2x2  system. 
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shows  that  the  relative  dynamic  gain  array  elements  are  the  limiting  values  of 

the  ratio  of  direct  transmittance  to  the  sum  of  direct  and  parallel  transmittance 
for  a  system  with  one  loop  open  and  the  rest  closed  Therefore  the  relative 

dynamic  gains  also  only  indicate  whether  problems  are  likely  to  be  experienced 
in  design  of  a  multiloop  control  system  for  a  particular  plant,  but  as  the 
examples  in  section  2.5  will  demonstrate  no  indication  of  closed  loop  interaction 

is  given  by  the  RDGA.  A  relationship  similar  to  equation  2.15  exist  for  higher 

order  systems 
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The  extended  interaction  quotient  and  the  diagonal  elements  of  the 
relative  dynamic  gain  array  are  related  as  follows 


1 


or 


Kx  = 


Mi/ 


(2.16) 


since 


1 1 


1-  X; 


■  — 

Hu 


m 


m 


n 


=  qr,*q  =  A1-;; 


j=i 

j*i 


Vi  % 


Cijq-y 


detQ 


(2.17) 


where  Cj.(s)  are  the  (i,j)'th  minor  or  Q(s).  The  information  contained  in  the  RDGA 
is  best  seen  by  plotting  it  as  an  array  of  polar  plots  similar  to  the  DNA 
display.  From  the  above  relationship  with  the  extended  interaction  quotient,  and 
the  interpretation  of  the  interaction  quotient,  it  is  evident,  that  the  RDGA  display 
should  be  interpreted  as  follows:  If  the  magnitude  of  the  diagonal  elements  are 
close  to  unity,  and  the  magnitude  of  all  other  elements  are  small  it  will  be 
possible  to  design  a  multiloop  control  system  to  give  satisfactory  control. 

Pairing  of  variables  with  magnitudes  much  different  from  unity  should  be 
avoided,  and  generally  elements  with  negative  real  part  should  not  be  placed  on 

the  diagonal,  since  the  corresponding  control  loop  will  show  non-minimum  phase 

properties. 


2.4.3  Tung  and  Edgar's  approach. 

Tung  and  Edgar  (1977,1978)  suggest  basing  the  variable  pairing  upon 
open  loop  step  responses  with  a  compensator  K^s)  =  G(0)  V  The  controlled 
and  manipulated  variables  are  to  be  paired,  so  the  dominant  coefficients  of  the 
time  response  arise  from  the  diagonal  elements  of  the  transfer  function  matrix 
G(s)G(0P1  However,  this  procedure  means,  that  in  systems  without  high 
frequency  interaction  the  variables  will  be  paired  so  the  first  output  is 
controlled  by  the  first  input  and  so  on,  i.e.  no  reordering  of  inputs  or  outputs 
The  resultant  pairing  can  be  considerably  different  from  a  similar  analysis 
performed  on  the  uncompensated  plant.  Tung  and  Edgar  do  not  mention  actually 
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implementing  the  precompensator  on  the  physical  system,  hence  the  final  con¬ 
troller  performance  could  be  unsatisfactory  for  systems  with  large  low 
frequency  interactions.  This  makes  Tung  and  Edgar's  approach  less  desirable  than 
the  RDGA. 


2.4.4  The  average  dynamic  gain  array. 

Gagnepain  and  Seborg  (1979)  also  suggest  using  open  loop  step 
responses  to  pair  variables.  They  extend  the  dynamic  potential  idea  of  Witcher 
and  McAvoy  (1977)  using  integrals  of  open  loop  step  responses  to  define  an 
average  dynamic  gain  array.  The  average  dynamic  gain  array  (ADGA)  is  defined 
by 


M(t)  =  D(t)o[D(t) 


(2.18) 


where  each  element  of  D(t)  is  an  integral  of  an  open  loop  step  response  from 
time  t1  to  time  t^  divided  by  the  integration  interval.  t1  is  choosen  as  the 
smallest  time  at  which  the  matrix  D(t)  is  not  singular  if  integration  was  from 
time  zero  to  time  tv  and  ta  is  choosen  so  the  integration  interval  t  =  t^  - 
1 1  is  equal  to  the  largest  time  constant  of  the  plant  transfer  function  matrix  . 
The  definition  of  the  ADGA  is  in  complete  analogy  with  the  definition  of  the 
RDGA,  and  it  therefore  has  many  of  the  properties  of  the  RDGA.  In  particular 
it  provides  an  answer  to  the  same  question  as  the  relative  gain  array,  i.e.  since 
D(t)  is  in  effect  a  matrix  of  average  open  loop  gains,  the  ADGA  measures  the 
significance  of  parallel  transmittance  in  the  same  way  and  for  the  same  situation 
as  the  RDGA.  The  difference  between  the  RDGA  and  the  ADGA  is,  that  the 
latter  tries  to  condense  the  dynamic  information  into  a  single  constant  matrix. 
Gagnepain  and  Seborg  show,  that  the  ADGA,  even  though  it  has  a  higher 
success  rate  than  the  steady  state  RGA,  fails  in  the  crucial  situations,  where 
stability  considerations  determine  the  best  pairing  of  variables. 


2.4.5  The  relative  transient  response  functions. 

Jaaksoo  (1979)  presented  an  exact  time  domain  equivalent  of  the  relative 
dynamic  gain  array.  Based  on  a  discrete  state  space  model  of  the  form 
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x(k+ 1 )  =  Ax(k)  +  Bu(k) 

(2.19) 

y(k)  =  Cx(k) 

relative  transient  response  functions  <pjj(k)  are  calculated.  These  functions  are 

defined  by 

k  .  .  k  .  , 

s:  C|A  bj  c IT. 

t=1  J  '  t=  1  J 

=“■ - =---- -  (2'20> 

Z  c.[Z 

where  Cj  is  the  ith  row  of  C,  b-  is  the  j'th  column  of  B  and  M  = 

l+Bj(CjBj)  1  Cj  in  which  Bj  is  B  without  the  j'th  column  and  Cj  is  C  without  the 

i’th  row.  Due  to  the  exact  equivalence  with  Bristol's  (1966)  approach  the  relative 
transient  response  functions  gives  a  measure  of  the  significance  of  parallel 
transmittance  with  one  loop  open  and  the  rest  closed,  and  not  an  indication  of 

closed  loop  interaction.  Jaaksoo  (1979)  suggests  using  only  the  arrays  obtained 
with  k=1  and  k  in  the  analysis.  This  corresponds  to  using  the  RDGA  with 

only  s  **♦  oo  and  s=0  respectively.  Jaaksoo  introduces  this  limitation  because  of 
difficulties  in  the  analytical  investigation  of  a  set  of  time  functions.  For  example 

the  matrix  power  series  in  equation  2.20  will  only  converge  if  all  eigenvalues 

of  A  and  (MjjA)  have  absolute  values  less  than  one,  and  the  denominators  can 

only  be  calculated  from  the  numerators  for  k=1,  cf.  relative  dynamic  gain  array. 
Further  drawbacks  of  the  relative  transient  response  functions  are:  it  is  not 
clear  how  the  function  should  be  displayed  or  interpreted  as  functions  of  the 

time  parameter  k;  and  they  require  a  significant  amount  of  numerical  calculations 

to  evaluate. 

2.4.6  Other  measures  of  interaction. 

None  of  the  measures  or  indices  discussed  so  far  require  any  knowlegde 
of  the  control  system  structure,  the  type  of  final  controllers  or  the  associated 
gains.  The  result  is,  that  the  indices  so  far  reviewed  are  more  of  warnings 

about  potential  difficulties,  than  indicators  of  closed  loop  interaction.  Three 
measures,  which  requirer  the  control  system  structure  and  the  final  controllers 
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to  be  known,  will  now  be  reviewed 

Davison  (1969)  suggests  a  non-minimum  phase  index  and  interaction  index 
based  on  a  state  space  model.  The  non-minimum  phase  index  is  evaluated  from 
eigenvalues  of  two  matrices  for  each  input  -  output  pair,  so  the  variable 

pairing  must  be  known,  and  the  only  type  of  controller  considered  is 

proportional  feedback.  The  interaction  index  is  the  difference  between  the 
closed  loop  non-minimum  phase  index  and  the  open  loop  non-minimum  phase 
index.  The  interaction  index  is  interpreted  in  terms  of  favourable  and 

unfavourable  interactions.  Since  as  has  been  shown  in  section  2.2  interaction 

transmittances  cannot  be  classified  as  good  or  bad  for  all  possible  load  or 

setpoint  changes  Davison’s  interpretation  is  wrong  The  interaction  index  can  only 

be  seen  as  a  measure  of  parallel  transmittance.  The  non-minimum  phase  index  is 

really  a  measure  of  the  importance  of  the  process  deadtime  relative  to  the 

dominant  time  constant,  and  as  such  has  little  to  do  with  interaction.  Both  the 
non-minimum  phase  index  and  the  interaction  index  are  of  limited  value  during 
control  system  design,  partly  due  to  the  calcuiational  effort  involved  in 
evaluating  them. 

Davison  and  Man  (1970)  suggest  an  interaction  index  based  on  the 

difference  between  a  single  closed  loop  response  and  the  multiloop  closed  loop 
response.  The  index  is  based  on  a  state  space  model  and  is  calculated  by 

solving  two  matrix  equations  iteratively  and  calculating  the  maximum  eigenvalues 
of  the  resulting  matrices,  for  each  closed  loop.  If  the  value  of  the  index  is 

small  for  a  given  control  loop  interaction  is  not  severe  in  that  loop  Even 
though  the  index  actually  measures  interaction  in  the  closed  loop  system  the 

amount  of  numerical  calculation  involved  in  evaluating  it  seems  prohibitive  for 

use  in  control  system  design.  Also  since  a  state  space  model  is  required  and 

the  complete  control  system  must  be  known,  it  appears  to  be  simpler  to 

conduct  some  closed  loop  simulations  in  order  to  evaluate  the  resultant  control 

system  performance  before  implementation. 

Suchanti  and  Fournier  (1973)  suggest  the  interaction  coefficients  defined 


for  each  loop  as 
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I.  =  E  (IE)  -  (IE)jl/(!E)j  (2.21) 

where  (IE)  *  is  the  error  integral  of  the  j'th  output  for  a  unit  step  change  in  the 

j'th  input  with  loop  j  closed  and  all  other  loops  open.  Similarly  (IE):  is  the  error 

J 

integral  of  the  j'th  output  for  a  unit  step  change  in  all  inputs  with  all  loops 

closed.  All  integrals  are  evaluated  with  proportional  plus  integral  controllers  tuned 

using  single  loop  techniques.  There  is  little  doubt,  that  the  index  gives  a 
measure  of  closed  loop  interaction,  and  it  is  possible  to  evaluate  the  coefficient 
for  any  control  system.  However,  calculation  of  the  index  requires  m+ 1  closed 
loop  simulations  for  a  mxm  system,  and  that  amount  of  closed  loop  simulation 
should  reveal  any  interaction  problems  without  the  need  to  calculate  error 
integrals.  Furthermore,  a  trial  and  error  approach  must  be  adopted  to  use  the 
coefficients  during  control  system  design,  i.e.  pairing  of  variables. 


2.5  The  DNA  -  an  alternative  to  measures  of  interaction. 

In  this  section  the  type  of  results  obtained  from  the  use  of  measures 
of  interaction,  such  as  the  interaction  quotient  of  Rijnsdorp,  the  RGA,  the  RDGA 
and  the  ADGA,  are  compared  with  the  information  obtained  from  the  DNA 
without  calculating  these  measures. 

Example  1: 


A  simple  2x2  transfer  function  matrix  will  show,  that  interaction  transmit- 
tances  can  be  severe  even  when  parallel  transmittance  with  one  loop  open  is 
small.  Consider  the  plant: 


Gp(s) 


1 

1s+  1 

_ 1_ 

_2s+1 


0.05 

IOs+1 

_ 1_ 

1  s+  1 


(2.22) 


and  assume  Ka(s)  =  H(s)  =  I  and  GL(s)  =  0.  Then  Rijnsdorps  interaction  quotient 
and  the  (2,2)-element  of  Bristol's  RDGA  are  respectively 

1 +2s+ Is*  , 1+12s+20s* 

K  =  0.05 - -  /***=  * - *  (2-23) 

1+12s+20s  0.95+  1  1.9s+  19  95s 


K  (s)  and  the  RDGA  M(s)  are  plotted  as  functions  of  frequency  in  figures  2.2 
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and  2.3  respectively.  The  magnitude  of  K(s)  is  small  at  all  frequencies,  and 
the  magnitude  of  the  diagonal  elements  of  M(s)  are  close  to  one  at  all 
frequencies.  Thus  both  measures  indicate  there  should  be  no  difficulty  in 
designing  a  multiloop  control  system,  which  is  expected,  since  the  TFM  in 
equation  2.22  is  almost  triangular. 

The  closed  loop  equation  for  the  second  output  with  two  proportional 
controllers  is 


y*  =' 


detF 


1 


1 


0.05 


*  /^+k2,k1(  X 

s+1  (s+1) 


1 + 1 2s+20s 


>r*+ki 


2,  '  * 


2s+  1 


(2.24) 


It  is  evident,  that  direct  transmittance  and  interaction  transmittance  have  very 
similar  dynamics  and  magnitude  Further,  the  parallel  transmittance  is  helpful  and 
significant  especially  at  low  frequencies.  With  the  second  loop  open  the 
equation  for  the  second  output  becomes 


o.05(s+i; 


s+  1 


s+  1 


*rz  kx“ 


(2s+  1)(  1 0s+  1)(s+  1  +k1 


,ra,+ki 


(2s+1)(s+1+k, 


(2.25) 


which  shows  insignificant  parallel  transmittance  under  these  conditions,  but 

comparable  direct  and  interaction  transmittances. 

Rijnsdorp's  interaction  quotient  and  the  RDGA  suggest  that  a  multiloop 
control  system  can  be  designed  without  difficulty,  but  they  fail  to  give  a 
reliable  measure  of  the  closed  loop  interaction.  For  this  particular  plant  the  DNA 
of  Gp(s)  gives  the  same  information  for  multiloop  design,  since  Gp(s)  have  no 
right  half  plane  poles  or  zeros  and  is  almost  triangular.  In  addition  the  DNA  of 
Gp(s),  shown  in  figure  2  4,  indicates  there  is  significant  closed  loop  interaction 
in  the  second  loop,  but  almost  none  in  the  first,  which  is  expected,  since  the 
TFM  is  almost  lower  triangular. 

Example  2: 

Gagnepain  and  Seborg  (1979)  consider  the  following  2x2  system 


Gp(s)  =  exp(-s) 


2  1.5 

IOs+1  1  s+1 

1.5  2 

1s+1  1 0s+ 1 


(2.26) 
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Figure  2.2  Polar  plot  of  Rijnsdorps  interaction  quotient  for  the  TFM  given  in 
equation  2.22  The  arrow  indicates  increasing  frequency. 


■ 
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Figure  2.3  Polar  plots  of  the  elements  of  the  relative  dynamic  gain  array  for 
the  TFM  in  equation  2.22.  Axis  are  only  shown  on  diagonal 
elements  and  the  origin  of  off  diagonal  elements  are  indicated  by  a 
plus  sign. 
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Figure  24  The  direct  Nyquist  array  for  the  TFM  given  in  equation  2.22 
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which  has  the  following  RGA  and  ADGA  respectively 


2.29 

-1.29" 

""-0.42 

1.42“ 

_- 1.29 

2.29 

1.42 

-0.42 

(2.27) 


The  two  measures  recommend  different  pairings.  Gagnepain  and  Seborg  (1979) 
show,  the  pairing  recommended  by  the  ADGA  is  structurally  monotonic  unstable 
according  to  a  theorem  due  to  Niederlinski  (1971).  The  failure  of  the  ADGA  is 
due  to  the  averaging  integration  process.  The  polar  plot  of  the  RDGA  in  figure 
2.5  suggests,  that  at  steady  state  the  first  input  variable  should  be  paired  with 
the  first  output  variable,  but  at  higher  frequencies  with  the  second.  The  DNA  of 
Gp(s)  leads  to  the  same  conclusion  as  seen  from  figure  2.6.  The  DNA  plot  not 
only  supplies  the  same  information  for  pairing  as  the  RDGA  but  also  provides  a 
basis  for  designing  the  appropriate  controller,  eg  using  the  DNA  design 
procedure. 

Example  3: 

Ahson  and  Nicholson  (1976)  and  later  -Taiwo  (1978)  consider  the  design 
of  a  compensator/controller  for  a  turbo-alternator  model  with  two  inputs  and 
two  outputs.  Simulations  by  Ahson  and  Nicholson  show  a  multiloop  control 
system  gives  very  unsatisfactory  control.  However,  the  steady  state  relative  gain 
array  is 


1.0235  -0.0235" 

.-0.0235  1.0235 


(2.28) 


which  indicates  there  should  be  no  difficulties  in  designing  a  multiloop  control 
system.  For  this  system  the  magnitude  of  the  elements  of  the  RDGA  change 
significantly  with  frequency,  as  seen  from  figure  2.7,  indicating  different  variable 

pairings  at  low  frequency  and  at  high  frequency.  The  same  information  is  also 
evident  from  the  DNA-plot  in  figure  2.8,  which  also  shows,  that  the  difficulties 
are  due  to  the  small  size  of  the  (1, 1  (-element  of  the  plant  transfer  function 

matrix.  The  transfer  function  model  used  by  Taiwo  (1978)  is  slightly  different 

from  the  one  used  by  Ahson  and  Nicholson  (1976).  The  RDGA  and  DNA  plots 

for  the  model  used  by  Taiwo  are  shown  in  figures  2.9  and  2.10  respectively. 
While  the  DNA  displays  of  figures  2.8  and  2.10  are  very  similar  and  readily 
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Figure  2.5  The  relative  dynamic  gain  array  for  the  TFM  given  in  equation  2.26 
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Figure  2.6  The  direct  Nyquist  array  for  the  TFM  given  in  equation  2.26 
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Figure  2.7  The  relative  dynamic  gam  array  for  turboalternator  using  model  of 
Ahson  and  Nicholson  (1976). 
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Figure  2.8  The  direct  Nyquist  array  for  turboalternator  using  model  of  Ahson 
and  Nicholson  (1976). 
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Figure  2.9  The  relative  dynamic  gain  array  for  turboalternator  using  model  of 
Taiwo  (1978). 
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Figure  2.10  The  direct  Nyquist  array  for  turboalternator  using  model 
(1978). 


of  Taiwo 
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reveal  the  difference  between  the  two  models  is  in  the  (2, 1)-element,  the 
RDGA  displays  of  figures  2.7  and  2.9  appear  to  come  from  two  highly 

different  models  The  DNA  displays  of  figures  2.8  and  2.10  further  show,  a 
full,  dynamic  compensator  (s)  is  probably  required  for  effective  control.  The 
complete  design  of  such  a  compensator  is  not  relevant  here,  but  note  that  the 

classification  of  transmittances  in  figure  2.1  provides  a  convenient  means  of 
comparing  some  of  the  design  alternatives: 

i.  Introduction  of  a  diagonal  precompensator  to  multiply  the  first  column 

of  Gp(s)  by  a  constant,  N1(  which  increases  the  direct,  parallel  and 
ri  ya  interaction  transmittance  by  N*. 

ii.  Introduction  of  a  diagonal  postcompensator  to  multiply  the  first  row  of 

Cf>{s)  by  a  constant,  Na,  which  increases  the  direct,  parallel  and  r^— >  y1 
interaction  transmittance  by  Na. 

iii.  For  a  given  increase  in  direct  transmittance  it  is  probably  best  to 
combine  the  above  two  approaches,  which  increases  the  direct  transmit¬ 
tance  by  N1  *N£,  the  parallel  transmittance  by  N^N^,  but  the  r1  — >  ya 
interaction  by  only  N^,  and  the  rJ^-*y1  interaction  by  only  Nx. 

Figure  2.11  shows  the  DNA  resulting  from  using  N1  =2  and  NA  =  5.  The  in¬ 
teractions  could  be  further  reduced  through  additional  design  steps.  Also  the 
diagonal  postcompensator  can  be  implemented  as  part  of  the  controller  and 
hence  causes  no  practical  difficulties.  Thus  it  is  seen,  that  the  DNA  provides  a 
better  basis  for  analysis  and  design  of  a  compensator  than  the  RDGA. 


2.6  Conclusions. 

Transmittances  in  a  closed  loop  multivariable  control  system  can 
advantageously  by  classified  as:  direct,  parallel,  interaction  and  disturbance  trans¬ 
mittances  Interaction  transmittances  in  an  mxm  system  can  further  be  expressed 
as  a  sum  of  1st,  2nd,  and  up  to  m-1'th  order  interaction  terms. 

Several  published  measures  of  interaction  were  reviewed.  The  interaction 
quotient,  the  RGA,  the  RDGA,  the  ADGA  and  the  relative  transient  response 
functions  have  the  following  shortcomings:  i)  they  are  based  on  a  system  with 


40 


Figure  2.11  The  DNA  for  turboalternator  using  model  of  Ahson  and  Nicholson 
(1976),  and  diagonal  pre-  and  postcompensators. 
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one  loop  open  rather  than  a  fully  closed  loop  system,  and  ii)  they  measure  the 
parallel  transmittance  with  all  but  one  loop  closed  rather  than  the  true  interac¬ 
tion  transmittance  They  are  more  of  a  warning  about  potential  difficulties  in 
designing  a  multiloop  control  system,  than  a  direct  quantitative  measure  of 
system  interaction. 

The  DNA  display  gives  an  indication  of  closed  loop  interaction  and  is  a 
suitable  tool  for  pairing  variables  as  the  first  step  in  a  multiloop  control  system 
design.  It  is  recommended  over  the  other  reviewed  measures. 


* 


3.  Matrix  dominance  and  transfer  of  dominance. 


3.1  Introduction. 

Multivariable  freauency  domain  design  techniques,  such  as  the  inverse 
Nyquist  array  (INA)  method  (Rosenbrock  (1969))  and  the  direct  Nyquist  array 
(DNA)  method  (Kuon  (1975)),  base  the  final  selection  of  controller  constants  and 
the  stability  analysis  on  only  the  diagonal  elements  of  respectively  the  inverse 
open  loop  transfer  function  matrix  (TFM)  and  the  open  loop  TFM.  In  order  to 
assess  stability  based  only  on  the  diagonal  elements  of  a  TFM,  some  condition 

must  be  put  on  this  TFM.  As  first  proposed  the  INA  method  required  the  open 

and  closed  loop  inverse  TFM's  to  be  diagonally  dominant,  and  the  DNA  method 
required  the  return  difference  matrix  to  be  diagonally  dominant,  in  the  sense 
defined  by  Rosenbrock  (1969).  The  application  of  the  multivariable  Nyquist  array 
(MNA)  techniques  to  industrial  problems  has  been  somewhat  limited  because  the 
diagonal  dominance  requirement  is  too  strict  for  many  industrial  systems,  where 
integrity  is  a  major  concern.  Recently  several  authors,  eg  Mee  (1976),  have 
explored  ways  of  relaxing  the  condition  for  application  of  the  multivariable 
Nyquist  array  stability  theorems 

In  this  chapter  past  work  on  extending  the  applicability  of  the  •  MNA 
stability  theorems  is  reviewed.  The  results  of  previous  work  on  M-matrices  are 
formalized  by  the  introduction  of  the  concept  of  matrix  dominance,  and  a 
necessary  and  sufficient  test  for  a  TFM  to  be  matrix  dominant  is  given.  A 

graphical  test  for  matrix  dominance  is  developed,  and  the  use  of  this  test  to 
aid  in  the  design  of  a  compensator  for  systems,  which  are  not  matrix 
dominant,  is  discussed 

The  implications  of  matrix  dominance  for  the  MNA  design  techniques  are 
demonstrated  Further,  it  is  shown,  that  a  row  dominant  system  can  be  made 

column  dominant  by  a  diagonal  similarity  transformation.  This  duality  is  finally 
exploited  in  a  selection  procedure  for  final  loop  gains  of  matrix  dominant 
systems. 
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3.2  Definition  of  and  test  for  matrix  dominance. 

In  the  extension  of  the  MNA  stability  theorems  the  concept  of 

M-matrices  plays  a  central  role.  The  next  subsection  contains  a  very  brief 

definition  of  an  M-matrix  and  some  associated  mathematical  tools.  This  is 
followed  by  sections  reviewing  previous  work  and  developing  tests  for  matrix 

dominance 

3.2.1  Mathematical  background  and  notation. 

A  real  m  x  m  matrix  A  =  { a 3  's  an  M-matrix  if  and  only  if  the 

VJ 

offdiagonal  elements  are  nonpositive  and  all  the  principal  minors  are  positive 

The  properties  of  M-matrices  has  been  extensively  treated  by  Ostrowski 

(1937,1956)  and  by  Fiedler  and  Ptak  (1962,  1966,  1967),  who  also  give 

alternative  criteria  for  judging  whether  or  not  a  real  square  matrix  is  an 

M-matrix.  (Fiedler  and  Ptak  use  the  term  matrix  of  class  K'  to  designate  an 

'M-matrix').  For  a  complex  m  x  m  matrix  Q  =  {qT$  the  companion  matrix 

(Begieitmatrix  (Ostrowski  (1956)))  B  =  £bij}  is  defined  as 

bj,-  =  lqi{l  and  b-,j  =  -lq-^1  i^j  (3.1) 

Further,  C  =  £c-.3  is  a  real  matrix  defined  as 

Cfi  =  0  and  c  —  =  Iqj^l  i^j  (3.2) 

and  W  =  diag[w;l  is  a  real  m  x  m  diagonal  matrix  with  w-  >  0. 

3.2.2  Previous  work  on  extending  MNA  stability  theorems. 

The  first  published  work  on  extending  the  class  of  TFM  to  which  the 

MNA  stability  theorems  are  applicable  is  due  to  Araki  and  Nwokah  (1975).  Araki 

and  Nwokah  use  the  concept  of  M-matrices  in  establishing  bounds  on  the 
transfer  function  h-(s)  between  the  i'th  input  and  the  i'th  output  of  a 

multivariable  system  with  the  i'th  loop  open  and  all  other  loop  closed  Their 

main  result  is,  that  for  a  open  loop  TFM  Q(s) 

lh.(s)  -  qjj(s)l  <  X(s)w;  (3.3) 

where  \{s)  is  the  largest  eigenvalue  of  CW”1  with  C  and  W  as  defined  in 

section  3.2.1.  The  authors  also  point  out,  that  any  one  bound  can  be  made 
arbitrarily  small  at  the  expense  of  making  another  bound  arbitrarily  large.  Finally 


< 


. 


44 


Araki  and  Nwokah  state  without  proof,  that  closed  loop  stability  of  an  open 
loop  stable  system  Q(s)  is  guaranteed  if  the  Nyquist  loci  of  q^s)  do  not 

encircle  the  critical  point  and  the  companion  matrix  of  the  return  difference 
matrix  F(s)  =  !+Q(s)  is  an  M-matrix.  Araki  and  Nwokah  do  not  clarify  the 

implications  of  their  findings  on  the  MNA  design  procedures,  and  give  no  simple 
way  of  ascertaining  whether  F(s)  has  a  companion  matrix,  which  is  an  M-matrix. 

Mee  (1976)  points  out,  that  a  practical  disadvantage  of  the  INA  technique 
is,  that  high  frequency  response  is  usually  not  well  known,  giving  rise  to 
difficulties  in  establishing  dominance  at  high  frequencies.  This  practical  difficulty 
is  also  reported  by  Leininger  (1979a)  and  Clement  (1980).  Mee  (1976)  shows, 
that  it  is  possible  to  test  stability  by  considering  the  Nyquist  loci  over  a  finite 
frequency  range  uo0  >  oo  >  -  instead  of  >  -  «=*=>  Mee 

also  introduces  the  use  of  diagonal  similarity  transformations  D(s)  =  diagfdj(s)}, 
d;(s)  real  and  positive,  to  make  a  nondiagonally  dominant  system  diagonally 
dominant  during  the  design  phase.  Complex  matrices  with  the  property  of  being 
diagonally  dominant  after  a  diagonal  similarity  transformation  are  matrices,  whose 
companion  matrix  is  an  M-matrix.  The  diagonal  similarity  transformation  shares 
dominance  among  or  transfers  dominance  to  different  rows  or  columns  The 

concept  of  dominance  sharing  is  also  discussed  by  Leininger  (1978,1979b),  who 
remarks,  that  applying  dominance  sharing  becomes  analytically  more  difficult  as 
the  order  of  the  transfer  function  matrix  increases.  Leininger  (1979b)  also 

restates  the  finite  frequency  range  result  of  Mee  (1976),  and  Leininger  (1979a) 
develops  a  function  minimization  procedure,  similar  to  pseudodiagonalization,  to 
design  a  constant  compensator,  to  achieve  diagonal  dominance.  The  dominance 
sharing  results  of  Mee  (1976)  and  Leininger  (1978,1979b)  are  essentially 

practical  applications  of  the  results  previously  stated  by  Araki  and  Nwokah 
(1975). 

Kantor  and  Andres  (1979)  use  a  method  identical  to  Araki  and  Nwokah's 

(1975)  to  calculate  normalized  Gershgorin  bands.  Kantor  and  Andres  calculate  the 
normalized  Gershgorin  radius,  ©  (s),  as  the  maximum  eigenvalue  of  CW”1  or 

W”1  C  with  w;  =  lqi(-(s)l,  where 


« 
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m 


(34) 


and  calculate  individual  radii  as  ©^(sllq^s)!.  The  djis)  ,  i 


m,  are  positive 


real  functions  of  s.  These  functions  are  the  elements  of  a  diagonal  similarity 
transformation  matrix  D(s).  Kantor  and  Andres  avoid  calculating  this  transformation 
by  transforming  the  problem  of  finding  minimal  Gershgorin  radii  into  an 
eigenvalue  problem.  The  largest  eigenvalue  must  be  choosen,  since  it  is  not 
possible  to  associate  eigenvalues  with  rows  or  columns  or  the  TFM  Q(s).  The 
main  result  of  Kantor  and  Andres  is,  that  Q(s)  is  similar  to  a  diagonally 

dominant  matrix,  i.e.  the  companion  matrix  of  Q(s'  is  an  M-matrix,  if  and  only 
if  ©’-(s)  ^  1  This  makes  it  possible  to  determine  from  a  simple  graph  of 
0^  (s)  versus  frequency  whether  or  not  the  companion  matrix  of  a  given  TFM 
is  an  M-matrix.  However,  if  the  test  is  negative  the  graphical  display  of  ©^(s) 
gives  no  guidelines  for  designing  a  suitable  compensator. 

Nwokah  (1980a)  formally  states  the  MNA  stability  theorems  for  TFMs, 
whose  companion  matrices  are  M-matrices.  Nwokah  (1980b)  points  out,  that  the 
results  of  Mee  (1976)  and  Kantor  and  Andres  (1979)  are  simple  consequences 
of  the  properties  of  M-matrices  Nwokah  (1980b)  also  misleadingly  labels 
complex  matrices,  whose  companion  matrix  is  an  M-matrix  as  Hadamard 

matrices.  The  term  'Hadamard  matrix'  is  generally  associated  with  matrices  with 
all  elements  -1  or  +1. 

3.2.3  Matrix  dominance. 

Since  the  terms  companion  matrix  and  M-matrix  do  not  clearly  show  the 
extension  of  the  diagonal  dominance  idea  of  Rosenbrockt  1 969),  and  the  term 
companion  matrix  in  control  is  often  defined  differently  from  Ostrowski’s  (1956) 

definition,  which  is  used  here,  it  is  proposed  to  call  the  extended  class  of 

systems  to  which  the  MNA  stability  theorems  are  applicable  matrix  dominant 
systems  Matrix  dominance  is  defined  as  follows 


Definition:  A  TFM  Q(s)  is  said  to  be  matrix  dominant  if  and  only 


if 


* 
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m  m 


(3.5) 


for  each  pair  (i,j),i£j. 


The  test  contained  in  this  definition  is  easily  programmed  on  a  digital  computer 
and  requires  only  a  few  more  calculations  than  a  diagonal  dominance  test.  It 
also  retains  the  sums  of  offdiagonal  elements  from  diagonal  dominance  test.  The 
above  definition,  which  is  adapted  from  a  theorem  of  Fiedler  and  Ptak  (1962), 
implies  that  the  companion  matrix  of  Q(s)  must  be  an  M-matrix  for  Q(s)  to  be 
matrix  dominant.  It  is  easily  seen,  that  the  class  of  diagonally  dominant  systems 
is  contained  in  the  class  of  matrix  dominant  systems,  i.e.  any  diagonally  dominant 
matrix  Q(s)  satisfies  inequalities  3.5 

The  property  of  matrix  dominance  is  identical  to  Mee's  (1976)  concept 
of  similarity  to  a  diagonally  dominant  system.  The  term  matrix  dominance  is 
suggested  here,  since  it  shows  the  relationship  to  diagonal  domi  nance,  and 
indeed  a  2  x  2  system  is  matrix  dominant  if  the  single  inequality 


iq<i  <s)Hq*»(s)l  >  Iq*!  (s)llqJZr(s)l 


(3.6) 


involving  all  elements  of  the  transfer  function  matrix  is  satisfied. 

3.2.4  Graphical  test  for  matrix  dominance. 

One  of  the  major  advantages  of  the  MNA  design  techniques  is  the  use 
of  graphical  displays  to  guide  the  designer  and  to  test  for  diagonal  dominance. 
A  graphical  test  for  matrix  dominance,  which  would  also  give  guidance  as  to 
the  design  of  a  compensator  for  a  system,  which  was  not  matrix  dominant, 
would  be  very  desirable.  The  inequalities  of  the  definition  for  matrix  dominance 
do  not  immediately  lend  themselves  to  a  graphical  test.  However,  a  slight 
rearrangement  of  the  inequalities  3.5  gives 


Hence,  a  system  is  matrix  dominant  if  the  real  non-negative  numbers  N:-  are 


*  M 

Therefore  a  set  of  plots 


less  than  one  for  all  pairs  (i,j),  i£j . 
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of  Nj.  as  a  function  of  frequency,  would  show  whether  or  not  a  system  is 

J 

matrix  dominant  For  a  2  x  2  system  only  one  graph  is  involved,  for  a  3  x  3 

system  3  graphs,  and  for  a  general  m  x  m  system  the  set  of  plots  contains 

m(m-1)/2  graphs.  This  is  less  than  the  number  of  plots  in  the  corresponding 

Nyquist  array,  but  more  than  in  a  graphical  test  for  diagonal  dominance  for  m 
>  3. 

To  illustrate  the  use  of  the  graphical  test  for  matrix  dominance  as  an  aid 

in  compensator  design  for  non-matrix  dominant  systems  a  transfer  function 

model  of  the  pilot  plant  double  effect  evaporator  in  the  Department  of 

Chemical  Engineering  at  the  University  of  Alberta  is  considered.  The  transfer 

function  model  calculated  from  the  continuous  fifth  order  state  space  model 

given  by  Fisher  and  Seborg  (1976)  is  shown  in  table  3.1.  A  set  of  three  plots 

to  graphically  test  F(s)  =  l+G(s)  for  column  matrix  dominance  is  shown  in  figure 

3.1,  and  the  similar  test  for  row  matrix  dominance  is  shown  in  figure  3.2.  In 
£  c  c 

figure  3.1  N1a,  and  NA2)  are  plotted  as  functions  of  frequency  for 

C  C 

frequencies  in  the  range  from  0.0038  to  7.73  rad/min.  N13>  and  N^  are  zero 
at  all  frequencies,  and  N1X  >  1  for  frequencies  lower  than  0.08  rad/min., 

hence  the  return  difference  matrix  is  not  column  matrix  dominant  at  low 
frequencies.  Figure  3.2  shows,  that  the  return  difference  matrix  is  also  not  row 
matrix  dominant  at  low  frequencies.  Since  only  Ni2>  is  different  from  zero  in 

figure  3.1,  this  indicates,  that  the  off  diagonal  elements  of  columns  1  and  2 
must  be  reduced  to  achieve  matrix  dominance  of  the  return  difference  matrix  at 
low  frequencies.  Figure  3.2  indicates,  that  the  offdiagonal  elements  of  row  3 

should  be  reduced  to  attain  matrix  dominance.  The  direct  Nyquist  array  (DNA) 
display  of  the  TFM  G(s)  is  shown  in  figure  3.3  (Note:  axis  are  only  drawn  on 
diagonal  elements,  and  the  origin  of  offdiagonal  elements  are  indicated  by  a  '+'). 
From  figure  3.3  it  appears,  that  the  open  loop  TFM  will  become  lower 
triangular  by  the  following  column  operation 

column  2  =  column  2  +  0.73*column  1  (3.8  ) 

which  will  tend  to  cancel  the  influence  of  element  (1,2)  and  at  the  same  time 
increase  the  direct  transmittance  (cf.  chapter  2)  in  the  second  loop,  and  slightly 
decrease  the  influence  of  element  (3,2).  After  the  above  compensation  the 


. 
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Figure  3.1  Graphical  test  of  pilot  plant  evaporator  return  difference  matrix  F(s) 
=  l+G(s)  for  column  matrix  dominance  Frequency  in  the  range  from 
0.0038  to  7.73  rad/min.  is  the  ordinate,  and  the  abscissa  the 
dimensionless  numbers  defined  by  equation  3.7. 
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Figure  3.2  Graphical  test  of  pilot  plant  evaporator  return  difference  matrix  for 
row  matrix  dominance  Frequency  range  0.0038  <  co  <  7.73 

rad/min. 
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Figure  3.3  Direct  Nyquist  array  display  of  pilot  plant  evaporator  TFM  G(s). 
Frequency  range  0.0038  <  ao  <  7.73  rad/min 


Table  3.1  Continuous  transfer  function  model  G(s)  of  pilot  plant  evaporator.  G(s) 
=  N(s)/d(s),  where  N(s)  is  a  polynomial  matrix  and  d(s)  is  the 

characteristic  polynomial 


d(s)  - 


n 


n 


n 


n 


n 


n 


n 


n 


n 


11 

\Z 

13 

21 

ZZ 

21 

31 

32, 

33 


(S) 

(s) 

(s) 

(s) 

(s) 

(s) 

(s) 

(s) 

(s) 


0.00025s* +0.036s3  +  0.82s4  +  1.0ss 
=  0  0013s* +0.032s3 
=  -0.00027s*  -0.032s3  -0.04 1  ss 
=  0.0 

=  -0.0000  1s*~0.00  13s3 -0.028sw 
=  -0.0000 2s- 0.00 28s*  -0.063sa  -0.077s4 
=  0.0 

=  -  0. 0000  1s-0. 00  15  s* -0.032s3 
=  0.00002s+0. 0029s*  +  0.065s3  +0.079s4 
=  -0.0000  1  S-0.00  1 4s*  -0.03  1  s3  -0.038s4 
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return  difference  matrix  is  row  matrix  dominant  as  shown  in  figure  3.4,  but  the 
return  difference  matrix  is  still  not  diagonally  dominant. 


3.3  Implications  for  MNA  design  techniques. 

In  this  section  it  is  shown,  how  the  property  of  matrix  dominance 
extends  the  applicability  of  the  MNA  stability  theorems.  The  following  theorem  is 
the  basis  for  applying  the  MNA  stability  theorems  to  matrix  dominant  systems: 
Theorem  3.7:  Let  the  square  complex  matrix  Q  be  matrix  dominant. 

Then  there  exist  a  diagonal  matrix  D  with  positive  diagonal 

^  4 

elements,  such  that  D  QD  is  diagonally  dominant. 

A  proof  of  this  theorem  is  given  in  appendix  C 

Since  l(l+D“1  QD)I  =  ID'1  (l+Q)DI  =  ID"1  ll(l+Q)IIDI  =  l(l+Q)l  the  return  difference 
determinants  of  F  =  l+Q  and  F  =  l+D"1  QD  are  the  same.  The  diagonal 

elements  of  Q  and  D”1  QD  are  also  identical,  and  it  is  well  known,  that  the 

eigenvalues  are  preserved  by  a  diagonal  similarity  transformation.  Hence  it 
follows  from  theorem  3  1  and  existing  results  for  diagonally  dominant  systems: 
Theorem  3.2:  Let  Q(s)  be  an  open  loop  transfer  function  matrix 
for  a  system  with  m-inputs  and  m-outputs.  Let  the  return  differ¬ 
ence  matrix  F(s)  =  l+Q(s)  be  matrix  dominant  for  all  s  on  the 
closed  right  half  plane  Nyquist  contour  D,  and  let  the  diagonal 

elements  of  Q(s)  map  D  into  V7;  ,  i=  1 m.  which  respectively 

encircle  the  (-1,0)  point  m  times  clockwise.  Then  the  closed  loop 

system  is  asymptotically  stable  if  and  only  if 
m 

E  n;  =  -p0  (3.9) 

i=  1 

where  p0  is  the  number  of  open  loop  system  poles  in  the  closed 
right  half  plane 

The  proof  follows  immediately  from  theorem  3  1  and  published  theorems  for 
diagonally  dominant  systems. 

A  similar  stability  theorem  in  terms  of  matrix  dominant  Q(s)’”1  and 
l+QIs)”1  can  be  stated  for  the  INA  design  procedure  However,  Q(s)  or  F(s) 


• 
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Figure  3.4  Graphical  test  of  pilot  plant  evaporator  return  difference  matrix  for 
row  matrix  dominance  after  compensation  Frequency  range  0  0038 
<  oo  7.73  rad/min. 
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being  matrix  dominant  does  not  imply  Q(s)”1  or  l+Qls)^1  is  matrix  dominant. 

To  illustrate  the  application  of  theorem  3.2  consider  the  following  2x2 
transfer  function 


G(s) 


1 

s+  1 


3  6 

L2  4.1 


(3.10) 


Rs)  =  l+G(s!  satisfies  inequalities  3.5  for  all  s,  so  F(s)  is  matrix  dominant  for  all 

s.  This  is  also  evident  from  the  graphical  test  shown  in  figure  3.5,  since  N1Jt 
<  1  at  all  frequencies.  However,  neither  G(s)  nor  Rs)  =  l+G(s)  are  diagonally 

dominant,  as  shown  by  the  Nyquist  array  display  in  figure  3.6.  Nevertheless 
theorem  3.2  allows  us  to  count  the  number  of  encirclements  of  the  (-1,0) 

point  by  qjj(s),  i=  1 ,2,  when  mapping  the  closed  right  half  plane  Nyquist  contour. 

The  number  of  encirclements  of  the  (-1,0)  point  is  zero,  and  since  G(s)  has  no 

poles  or  zeros  in  the  closed  right  half  plane,  the  closed  loop  system  is 
asymptotically  stable  according  to  theorem  3.2. 

Figure  3.6  also  illustrates  the  trade  off  between  matrix  dominance  and 
diagonal  dominance.  The  final  loop  gains  must  be  selected  so  the  characteristic 

loci  encircle  the  (-1,0)  point  clockwise  the  appropriate  number  of  times,  e.g. 

zero  times  for  a  system  with  no  right  half  plane  poles  or  zeros.  The 
characteristic  loci  are  known  to  lie  inside  the  Gershgorin  circles,  and  the  bigger 
the  Gershgorin  circles  the  larger  the  uncertainty  about  the  location  of  the 

characteristic  loci.  Because  of  this  uncertainty  the  gain  k;  of  the  i'th  loop  must 

be  choosen  conservatively,  so  —  1  /kj  lies  to  the  left  of  the  band  swept  out  by 

the  Gershgorin  circles  of  the  i'th  diagonal  element  (if  the  system  has  no  right 

half  plane  poles  or  zeros).  For  the  system  considered  here,  this  means,  a  high 

gain  could  be  applied  to  loop  1,  but  a  very  low  gain  would  have  to  be 

choosen  for  the  second  loop,  i.e.  a  conservative  gain  would  probably  be 
selected  for  the  second  loop.  Alternatively  a  trial  and  error  procedure  could  be 

adopted  in  which  gains  are  initially  selected  based  on  qjj(s),  and  the  stability 

test,  i.e  the  matrix  dominance  test,  repeated.  For  a  diagonally  dominant  system 

the  DNA-display  would  have  smaller  Gershgorin  circles  and  hence  there  would 
be  less  uncertainty  in  the  location  of  the  characteristic  loci,  and  therefore  the 
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Figure  3.5  Graphical  test  of  return  difference  matrix  of  system  in  equation  3.10 
for  matrix  dominance  Frequency  range  0. 1  <  oo  <.  10.0. 


. 
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Figure  3.6  Direct  Nyquist  array  of  2  x  2  TFM  G(s)  with  Gershgorin  circles 
superimposed  on  diagonal  elements.  Frequency  range  0.1  <  60  < 
10.0. 
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probability  of  choosing  a  conservative  gain  in  any  loop  would  be  smaller.  In 
section  3.5  a  procedure  based  on  transfer  of  dominance  is  developed  to  aid  in 
the  selection  of  final  loop  gains  for  matrix  dominant  systems. 


3.4  Duality  of  row  and  column  dominance. 

The  first  theorem  of  the  previous  section  showed  how  a  matrix 
dominant  system  can  be  made  diagonally  dominant  by  a  diagonal  similarity  trans¬ 
formation.  In  this  section  it  is  shown,  that  the  diagonal  similarity  transformation 
can  be  designed  to  make  the  matrix  dominant  system  either  row  diagonally 

dominant  or  column  diagonally  dominant.  The  following  theorem  establishes  this, 
since  row  diagonal  dominance  and  column  diagonal  dominance  are  both  special 
cases  of  the  more  general  matrix  dominance. 

Theorem  3.3:  Row  and  column  diagonal  dominance  are  dual  in  the 

sense,  that  if  a  complex  square  matrix  Q  is  column  diagonally 
dominant,  then  there  exist  a  diagonal  matrix  D  with  positive  diago¬ 
nal  elements,  such  that  D‘"’*QD  is  row  diagonally  dominant,  and 

vise  versa. 

The  proof  is  given  in  appendix  C. 

The  results  of  theorems  3.1  and  3.3  immediately  lead  to  the  following 
corollary: 

Corollary  3.7:  Diagonal  similarity  transformations  to  make  a  matrix 
dominant  system  either  row  or  column  diagonally  dominant  always 

exist. 

The  similarity  transformation,  which  makes  a  matrix  dominant  system  row 

diagonally  dominant,  will  in  general  be  different  from  the  similarity  transforma¬ 
tion,  which  makes  the  same  matrix  dominant  system  column  diagonally  dominant. 


The 

duality 

of  row 

and 

column 

dominance  can 

be  viewed 

as 

further 

support  of 

the 

idea,  that 

the 

dominance  requirements 

of  the  MNA 

design 

techniques 

are 

artificial 

requirements 

posed  on  the 

system 

for 

purely 

mathematical  reasons.  In  general  dominance,  particularly  diagonal  dominance,  is  not 
a  necessary  requirement  for  a  satisfactory  control  system  design,  although  a 
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highly  dominant  system  often  means  a  system  with  little  interaction. 

The  use  of  diagonal  similarity  transformations  as  an  aide  in  controller 
design  by  transfering  dominance  is  considered  in  the  following  section. 


3.5  Transfer  of  dominance. 

The  final  proportional  loop  gains  of  a  matrix  dominant  multivariable 
system  could  be  selected  on  the  basis  of  the  direct  Nyquist  array  display  of 
such  a  system  with  its  associated  wide  Gershgorin  bands.  However,  such  a 


procedure 

could  as 

discussed  in 

section  3.3 

lead  to  the 

selection 

of 

conservative 

gains  or 

loose  control 

in  one  or 

more  loops. 

In  order 

to 

overcome  this  deficiency  a  diagonal  similarity  transformation  could  be  used  to 
make  the  matrix  dominant  system  diagonally  dominant  during  the  design.  This 
approach  was  first  brought  forward  by  Mee  (1976),  who  also  pointed  out,  that 
the  similarity  transformation  need  not  be  implemented  in  the  field,  because  it 
does  not  change  any  transmittances,  but  only  redistributes  dominance,  taking 
from  highly  dominant  rows  and/or  columns  and  giving  to  non-dominant  rows 
and/or  columns.  Mee  (1976)  considered  only  2x2  systems,  where  the 
redistribution  or  sharing  of  dominance  was  simple,  but  gave  no  procedure  for 
calculating  the  necessary  similarity  transformation  for  a  generai  m  x  m  system. 
An  algorithm  for  calculating  the  required  similarity  transformation  for  a  general 
m  x  m  system  will  be  outlined  below. 

The  similarity  transformation  of  theorem  3.3  can  be  found  as  the  solution 
to  inequalities  9.22  and/or  9.23  of  appendix  C.  The  solution  of  these 
inequalities  can  be  cast  into  a  linear  programming  problem.  Consider  the  m  x  m 
transfer  function  matrix  Q(s),  which  is  assumed  to  be  matrix  dominant.  The 
diagonal  similarity  transformation  to  maximize  the  row  diagonal  dominance  is 
found  by  maximizing  the  following  penalty  function 


J  =  dm>+1  (s)  (3.1  1 ) 

subject  to  the  constraints 

d;(s)  >0  i  —  1 . m+1  (3.12) 


and 
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m 


lqjj(s)ldj(s)  -  21  lq;j(s)ldj(s)- 


(s)dm-M(s)  >  0 


(3.13) 


for  i 


m,  and 


d1  (s)  <  1 


(3.14) 


where  /3 ;  (s)  are  row  weighting  factors.  These  weighting  factors  allow  the 
designer  to  make  certain  rows  strongly  dominant  at  the  expense  of  other  rows, 
i.e.  dominance  is  transfered  between  rows.  An  equivalent  algorithm  for  transfer 
between  columns  is  easily  formulated,  in  which  case  <3;(s)  are  column 
weighting  factors.  The  optimization  problem  outlined  above  can  be  seen  to  have 
the  structure  of  a  general  linear  programming  problem 

maximize  J  =  g(x)  (3.15) 


subject  to 


A  x  ^  0 


(3.16) 


which  is  easily  solved  using  available  standard  procedures,  such  as  the  simplex 
algorithm.  The  solution  is  done  numerically  at  a  set  of  discrete  frequencies  s  = 
ioo^  ,  k  =  1 . ,p. 

The  algorithm  allows  certain  rows  or  columns  to  be  made  strongly  row 

or  column  diagonally  dominant,  thus  decreasing  the  uncertainty  about  the  stability 
margin  for  certain  important  input/output  pairs.  This  might  be  an  important 

consideration  in  many  practical  process  control  applications.  The  above  algorithm 
also  removes  any  analytical  difficulties  in  sharing  dominance  as  the  order  of  the 

transfer  function  matrix  increases,  since  the  designer  simply  assigns  weighting 
factors  according  to  the  importance  of  each  output  (row)  or  input  (column). 

The  flexibility  of  the  outlined  algorithm  as  compared  to  the  method  of 
Kantor  and  Andres  (1979)  for  calculation  of  minimal  Gershgorin  band  radii  is 

demonstrated  by  considering  a  constant  2x2  matrix.  Consider  the  matrix 


f3  6' 


1  4 


(3.17) 
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Then  the  eigenvalues  of  CW~^  with  C  and  W  as  defined  in  section  3.2.1  are 
-+1/VZ  and  the  method  of  Kantor  and  Andres  (1979)  gives  the  following 
diagonally  dominant  matrix 


3  2.12 

2.83  4 


(3.18) 


The  linear  programming  scheme  developed  in  this  section  gives  different  results 
depending  on  the  choice  of  the  weighting  factors.  For  /34  =  /3^  =  1  the 

transformed  system  becomes 


3  2.50" 

2.40  4 

U 

or  for  /31  =  2  and  /2>a  =  1  the  transformed  system  becomes 


(3.19) 


3  2.14 

2.80  4 


(3.20) 


which  is  almost  identical  to  3.18.  The  higher  weighting  /3>1  gives  a  smaller 
Gershgorin  radius  for  the  first  row.  Hence  the  example  demonstrates,  that  with 
the  transfer  of  dominance  algorithm  developed  here  the  designer  has  the 
freedom  to  make  the  bands  on  certain  important  loops  small  at  the  expense  of 
larger  bands  in  other  loops,  i.e.  dominance  can  be  arbitrarily  transfered  from 
one  row  or  column  to  another  by  the  choice  of  appropriate  weighting  factors, 
<Z  \  (s). 


3.6  Conclusions. 

The  past  efforts  in  extending  the  class  of  transfer  function  matrices  to 
which  the  MNA  stability  theorems  are  applicable  were  reviewed,  and  the 
common  property  of  this  class  of  matrices  was  defined  as  matrix  domi  nance. 
A  graphical  test  for  matrix  dominance,  that  also  aides  in  compensator  design 
was  developed,  and  its  application  illustrated. 

The  implications  of  matrix  dominance  for  the  MNA  design  techniques 
were  discussed,  and  row  and  column  diagonal  dominance  shown  to  be  dual  with 
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respect  to  a  diagonal  similarity  transformation. 

A  systematic  algorithm  for  transfering  dominance  in  tranfer  function 
matrices  of  any  size  were  developed,  and  its  flexibility  illustrated  by  applying  it 
to  calculation  of  minimal  Gershgorin  radii.  The  transfer  of  dominance  algorithm 
developed  as  part  of  this  study  gives  the  designer  more  freedom,  than  other 
published  procedures  for  calculating  minimal  Gershgorin  radii. 


4.  Nyquist  exact  loci  (NEL)  design  technique. 


4.1  Introduction. 

In  the  past  ten  years  the  classical  frequency  domain,  single-input, 
single-output  control  system  design  techniques  of  Nyquist  and  Bode  have  been 
extended  to  the  design  of  multi-input,  multi-output  control  systems.  Many 

significant  contributions  to  the  development  of  frequency  domain  control  system 

design  techniques  have  been  collected  into  a  book  edited  by  MacFarlane  (1979). 
Kuon  (1975)  reviews  the  developments  in  multivariable  frequency  domain  design 

techniques.  The  major  differences  between  multivariable  frequency  domain  design 

techniques,  such  as  the  multivariable  Nyquist  array  methods,  the  characteristic 
loci  method,  and  classical  single  variable  techniques  are  discussed  in  the 
following  paragraphs. 

The  first  major  extension  of  frequency  domain  design  techniques  to 

multivariable  systems  was  the  development  of  the  inverse  Nyquist  array  (INA) 

design  procedure  by  Rosenbrock  (1969).  The  inverse  open  loop  and  inverse 
closed  loop  TFM  must  be  dominant  before  system  stability  can  be  evaluated. 
Dominance  can  (usually)  be  achieved  by  performing  row  or  column  operations  on 
the  set  of  Nyquist  plots  -corresponding  to  the  m  x  m  elements  of  the  open 

loop  TFM.  Once  the  desired  degree  of  dominance  has  been  achieved,  stability 
can  be  determined  by  counting  the  encirclements  of  the  origin  and  the  critical 

point  by  the  Nyquist  loci  of  the  m  diagonal  elements  of  the  open  loop  TFM 
and  the  corresponding  Gershgorin  or  Ostrowski  circles  Because  of  the  inverse 
nature  of  the  plots  and  the  uncertainty  introduced  by  the  Gershgorin  or 
Ostrowski  bands,  the  frequency  and/or  time  domain  relationship  between  a 
specific  pair  of  input  -  output  variables  is  not  as  clear  as  with  the  classical 

SISO  Nyquist  plots. 

The  second  major  extension  of  frequency  domain  control  system  design 

techniques  was  the  introduction  of  the  characteristic  loci  (CL)  design  procedure 

by  MacFarlane  and  Belletrutti  (1973).  The  CL  technique  uses  a  set  of  plots  of 
the  eigenvalues  of  the  open  loop  TFM  to  assess  stability,  and  uses  both  the 
eigenvalues  and  eigenvectors  of  the  open  loop  TFM  to  design  a  multivariable 
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compensator/controller.  MacFarlane  and  Bel  letrutti  avoid  putting  any  conditions, 
such  as  dominance,  on  the  open  loop  TFM  by  using  the  eigensystem  of  this 

matrix  in  the  design  procedure.  The  designer,  using  the  CL  technique,  only  has 
to  work  with  m  diagrams  for  a  m-input,  m-output  multivariable  system. 

However,  these  m-diagrams  have  no  one  to  one  relationship  with  the  simple 
input/output  relationships  of  SISO  systems,  since  eigenvalues  cannot  in  general 
be  associated  with  specific  input/output  pairs.  Furthermore  the  familiar  SISO 

control  system  design  terms,  such  as  gain  margin,  phase  margin,  bandwidth,  rise 
time,  overshoot  and  settling  time  change  meaning  or  cannot  be  evaluated  directly 

from  the  characteristic  loci  diagrams. 

Other  significant  extensions  of  classical  frequency  domain  techniques  to 
the  design  of  multivariable  control  systems  include  the  sequential  return 

difference  (SRD)  design  procedure  proposed  by  Mayne  (1972,  1973,  1974),  and 
the  direct  Nyquist  array  (DNA)  design  procedure  proposed  by  Rosenbrock  (1974) 

and  described  in  detail  by  Kuon  (1975).  The  SRD  technique,  as  the  name 
suggests,  constructs  the  multivariable  controller  from  a  sequence  of  single  loop 
designs.  The  i'th  stage  of  this  sequence  designs  a  controller  for  the  i'th  loop 
taking  into  consideration  the  already  completed  designs  of  loops  1  to  i-1.  The 
drawback  of  the  SRD  procedure  is,  that  the  final  design  depends  on  which 

sequence  the  loops  are  closed  in,  and  no  general  method  exists  for  choosing 

the  best  sequence  The  DNA  technique  is  similar  to  the  INA  technique,  but  uses 
the  open  loop  TFM  directly  to  design  a  multivariable  compensator/controller,  ana 
requires  the  return  difference  matrix  to  be  dominant  before  system  stability  can 
be  evaluated  If  the  return  difference  matrix  is  strongly  diagonally  dominant 

familiar  SISO  control  system  design  specifications  can  be  evaluated  directly  from 

the  diagonal  elements  of  the  direct  Nyquist  array.  However,  in  the  general  case 

the  meaning  of  terms,  such  as  rise  time,  overshoot  and  settling  time,  still 
remain  fuzzy,  since  the  diagonal  elements  of  the  DNA  do  not  represent  the 
complete  transmittance  between  input  -  output  pairs  with  all  other  loops  closed, 
but  only  the  direct  part  (cf.  chapter  2). 

In  this  chapter  a  practical  and  systematic  multivariable  frequency  domain 
design  procedure,  called  the  Nyquist  exact  loci  (NEL)  technique,  is  outlined.  The 
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NEL  design  procedure,  which  is  applicable  to  a  large  class  of  industrial  systems, 

uses  the  transfer  functions  hj(s),  i  =  1 ,m  cf.  figure  1.2,  which  represents  the 

complete  transmittance,  direct  and  parallel,  between  the  i'th  input  and  the  i'th 
output  with  the  i'th  loop  open  and  all  other  loops  closed.  The  transfer 
functions  h;(s)  are  direct  multivariable  parallels  of  the  classical  SISO  transfer 
function  g(s).  The  transfer  functions  h5(s)  are  first  introduced  in  the  next  section 
and  the  advantages  of  working  with  them  discussed.  This  is  followed  by  a 
general  step  by  step  outline  of  the  NEL  design  technique.  The  different  steps 
of  the  NEL  procedure  are  then  discussed  in  detail  from  a  theoretical  and  a 
practical  viewpoint.  Finally  the  simultaneous  gain  calculation  algorithms,  which  are 
a  central  part  of  NEL,  are  described. 


4.2  The  Nyquist  exact  loci  of  h;{s). 

The  central  feature  of  the  NEL  procedure  is  the  use  of  the  transfer 
functions  hj(s)  ,  i  =  1 m,  to  simultaneously  calculate  the  constants  of  m  single 
loop  controllers  and  to  determine  closed  loop  stability.  The  transfer  function 
hj(s)  represents  the  complete,  direct  and  parallel,  transmittance  between  the  i’th 
input  and  the  i'th  output  when  the  i'th  loop  is  open  and  all  other  loops  are 
closed.  It  is  thus  the  exact  single-input,  single-output  system  for  which  the  i'th 
controller  should  be  designed  to  account  for  the  multivariable  nature  of  the 

system.  The  i'th  controller  is  the  i'th  diagonal  element  of  K^(s),  cf.  figure  1.1, 

and  the  transfer  function  h;(s)  is  the  multivariable  analogue  of  the  S'SO  transfer 
function  g(s),  cf.  figure  1.2.  It  is  therefore  clear,  that  the  transfer  function  h-((s) 
has  physical  meaning,  and  hence  terms  familiar  from  classical  single  variable 

design,  such  as  gain  margin,  phase  margin,  bandwidth,  rise  time,  overshoot  and 

settling  time,  can  be  evaluated  directly  and  meaningfully  from  Nyquist  plots  of 
h-(s). 

The  transfer  functions  hj(s)  are  nonlinear  functions  of  the  elements  of  the 
compensated  plant  TFM  Q(s)  and  the  elements  of  the  diagonal  controller  matrix 
(s).  The  following  general  expression  for  h.(s)  is  derived  in  appendix  A 

5v  I 
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h-,(s) 


qii(s) 


m 

e: 

j=i 

j^i 


<f)jj  (s)kj(s)q.j(s)/kj(s)  = 


q;:(s) 


m 

r 

j=i 


<t>ji  |s>qji(s) 


(4  1 


where 

a. 

cofactor 

of 

element  (i,j)  of 

F(s) 

(4.2) 

<Pij(s) 

cofactor 

of 

element  (i,i)  of 

F(s) 

and  F(s)  = 

l+Q(s)Ka(s) 

= 

l+G(s)K4  (s)K^(s). 

The  above  expression  for  i 

=  m  gives 

the  same  transfer  function  as  an  expression  given  by  Rosenbrock  (1972). 
However,  for  i  =  m  the  expression  given  by  Rosenbrock  is  the  transfer 

function  with  loops  1 . ,i—  1  closed  and  loops  i . m  open  Such  a  function  has 

been  used  by  Mayne  (1974)  in  his  SRD  design  procedure.  The  advantage  of  the 
above  expression  for  h*(s)  is,  that  it  facilitates  simultaneous  design  of  the  m 
controllers  in  K^(s)  =  diag{kj(s)}. 

Equation  4.1  results  in  the  following  expression  for  h^(s)  for  a  2  x  2 
system  with  loop  1  closed  (cf.  figure  1.1  and  1.2): 

k  i<s)q2/1  (s) 

ha(s)  =  q^(s)  + - q1Jt(s)  (4.3) 

1+k1(s)q11  (s) 

and  similarly  h^s)  with  loop  2  closed  is  given  by: 

ka (s)q1a  (s) 

hj(s)  =  qn(s)  + - qzi(s)  (4.4) 

1+k^(s)q^(s) 

Since  the  transfer  function  h.(s)  is  the  exact  transmittance  for  which  k«(s)  should 
be  designed  to  take  into  consideration  the  multivariable  nature  of  the  plant,  the 
polar  plot  of  h;(s)  for  s  on  the  Nyquist  £)-contour  is  called  a  Nyquist  exact 
locus,  and  the  design  procedure  based  on  these  loci  is  called  the  Nyquist  exact 
loci  (NEL)  technique.  The  NEL  design  procedure  is  outlined  in  the  next  section. 


4.3  The  Nyquist  exact  loci  (NEL)  design  procedure. 

The  NEL  design  procedure  for  the  design  of  a  multivariable  controller  for 
a  general  multivariable  plant  involves  the  following  steps,  which  are  discussed  in 
greater  detail  in  the  following  subsections: 
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i.  Obtain  a  transfer  function  model  description  of  the  plant  or 

system. 

ii.  Check  if  the  return  difference  matrix  is  matrix  dominant  with  (s) 

=  K^(s)  =  I. 

If  yes  then  continue  with  the  design  of  K^(s)  leaving  K1(s)  =  I,  or 

design  K  (s)  to  reduce  interaction, 
else  design  (s)  and/or  reduce  magnitude  of  elements  of 
K^(s)  to  make  F(s)  matrix  dominant. 

iii.  Choose  a  set  of  initial  gains,  kj,  i  =  1 m. 

iv.  Plot  the  Nyquist  exact  loci  of  h5(s),  i  = 

v.  Specify  stability  margins  for  each  loop  based  on  the  Nyquist  exact 
loci. 

vi.  Simultaneously  calculate  controller  constants,  which  satisfy  the 

stability  margins  of  all  loops,  cf.  section  4.4. 

vii.  Check  if  F(s)  is  matrix  dominant. 

viii.  Determine  if  the  design  is  satisfactory  by  plotting  the  Nyquist 

exact  loci. 

If  yes  and  F(s)  is  matrix  dominant,  then  the  design  is  completed 

after  checking  stability  using  the  Nyquist  exact  loci, 

else  if  yes  and  F(s)  is  not  matrix  dominant,  then  check 
stability  by  method  which  does  not  require  system 

dominance,  eg  characteristic  loci. 
else  relax  design  specification  and/or  improve  design  of 

compensator  K^(s)  to  meet  the  design  specifications. 

The  flow  chart  in  figure  4  1  summarizes  the  NEL  design  procedure,  and  the 

following  subsection  discusses  each  step  in  detail. 

4.3.1  Transfer  function  model  description. 

The  NEL  procedure  works  with  a  TFM  model  of  the  system  under 

consideration  as  the  basis  for  the  control  system  design.  The  TFM  model  can 

be  either  a  continuous  model  G(s),  a  discrete  model  G(z),  or  a  bilinearly 
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Figure  4.1  Flow  chart  of  NEL  design  procedure  showing  the  different  steps  and 
design  loops. 
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transformed  discrete  model  G(w).  If  a  system  description  is  available  in  some 
other  form,  e.g  a  partial  differential  equation  model  or  a  state  space  model,  a 

transfer  function  model  must  be  calculated  from  this  form  before  the  NEL 

procedure  can  be  used.  In  a  practical  computer  aided  implementation  of  the  NEL 
procedure  the  TFM  model  is  evaluated  at  a  set  of  discrete  frequency  points  on 
the  Nyquist  D~contour,  and  all  numerical  calculations  then  proceed  from  the 
resulting  set  of  complex  matrices,  rather  than  using  algebraic  manipulations  of 
polynomials  in  s,  z  or  w.  The  Nyquist  D-contour,  along  which  the  transfer 
functions  are  evaluated  depends  on  the  type  of  TFM  model.  A  discussion  of 
the  different  Nyquist  contours  in  the  s-plane,  z-plane  and  w-plane  is  given  in 

appendix  D.  Theoretically  the  transfer  functions  must  be  evaluated  along  the 

appropriate  complete  Nyquist  contour,  however  results  by  Mee  (1976)  and 

Leininger  (1979b)  show,  that  it  is  possible  to  use  only  a  finite  frequency  range 

0  <  co  <  co0  in  the  Nyquist  analysis,  provided  the  transfer  function  matrix 

is  bounded  for  all  frequencies  above  .  The  use  of  a  finite  frequency 

range  makes  the  use  of  the  NEL  design  procedure  and  other  multivariable 

Nyquist  array  techniques  more  practical,  since  it  is  considerably  easier  to  make 

e.g.  the  return  difference  matrix  dominant  over  a  finite  frequency  range,  than 
over  the  complete  Nyquist  £>-contour. 

Many  control  systems,  especially  if  K^s)  £  I,  are  implemented  using  a 

digital  computer.  Such  implementations  normally  involve  sampling  of  inputs  and 

outputs  at  equidistant  instants  of  time,  and  holding  the  plant  inputs  constant 

between  sampling  instants.  Hence  a  continuous  plant  is  controlled  by  a  discrete 
controller.  The  design  of  this  discrete  controller  can  be  approached  in  several 
ways: 


Based  on  a  continuous  TFM  G(s),  a  compensator/controller  K1  (s)K^(s) 
is  designed  in  the  s-plane.  The  discrete  implementation  is  then 

obtained  by  mapping  the  poles  and  zeros  of  K1(s)K^(s)  into  the 

z-plane  to  obtain  K^(z)K7t(z),  and  deriving  a  difference  equation 

from  K1(z)K2<(z)  for  time  domain  implementation.  This  procedure  is 
theoretically  incorrect,  and  as  demonstrated  by  Franklin  and  Powell 
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(1980)  for  SISO  system  it  does  not  give  a  control  system,  which 
is  in  agreement  with  the  design  specifications. 

ii.  From  the  continuous  TFM  G(s)  a  discrete  TFM  G(z)  is  calculated  by 

z-transformation  of  G(s)G^(s),  where  GjJs)  is  the  transfer  function 
of  a  zero  order  hold.  Then  a  compensator/controller  K^(z)Kx{z)  is 
designed  in  the  z-plane,  and  the  difference  equation  obtained 

directly.  This  procedure  is  theoretically  correct,  and  the  discrete 
controller  is  as  good  as  the  model  of  the  system.  However,  the 
z-plane  analysis  is  different  from  the  s-plane  analysis,  which  most 
designer  are  more  familiar  with.  Also  the  discrete  transfer 

functions  are  not  rational  functions  of  frequency. 

iii.  From  the  continuous  TFM  G(s)  a  bilinear ly  transformed  discrete  TFM 

G(w)  is  calculated  by  z-transformation  of  G(s)G^_(s)  to  give  G(z), 
followed  by  bilinear  transformation  of  G(z)  to  give  G(w).  Then  a 
compensator/controller  K^wlK^w)  is  designed  in  the  w-plane  and 

the  inverse  bilinear  transformation  used  to  obtain  K<(z)Ka(z),  from 
which  the  difference  equation  is  derived.  This  procedure  is 
theoretically  correct,  and  the  resulting  controller  is  as  good  as  the 
model  of  the  system.  The  advantage  of  performing  the  exact 
bilinear  transformation  from  the  z-plane  to  the  w-plane  is,  that  the 
analysis  in  the  w-plane  is  very  similar  to  the  analysis  in  the 
s-plane,  eg  stable  poles  and  zeros  lie  in  the  left  half  plane,  and 
the  transfer  functions  are  again  rational  functions  of  frequency,  cf. 
appendix  D 

iv.  Based  on  a  continuous  TFM  G(s)Gk(s)  a  compensator/controller 

K1(s)K2/(s)  is  designed  in  the  s-plane.  The  discrete  implementation  is 
then  obtained  by  mapping  the  poles  and  zeros  of  Ki  (s)K^(s)  into 
the  z-plane  to  get  K1(z)Ka(z)  from  which  a  difference  equation  is 
derived  for  implementation.  This  procedure  is  not  theoretically 
rigorous,  but  as  discussed  below,  if  the  discrete  transfer  functions 

are  calculated  from  the  continuous  using  the  Tustin  rule,  then  this 
procedure  is  identical  to  iii)  given  above  However,  readers  should 


. 
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be  aware  of  some  distortion  in  mapping  from  the  s  to  the 

z-plane  (Franklin  and  Powell  (1980))  in  using  this  approach.  In  the 
overall  sence  this  approach  is  a  good  practical  approximation, 
which  accounts  for  the  major  feature  of  the  discrete  control 

systems  implementation,  the  zero  order  hold.  An  alternative  to  this 
approach  is  the  hold  equivalence  method  also  described  by  Franklin 

and  Powell  (1980). 

The  four  approaches  to  design  of  a  discrete  controller  system  for  a  continuous 

plant  are  summarized  in  table  4.1.  If  a  discrete  TFM  G(z)  or  G(w)  is  available 

the  initial  transformations  of  approaches  ii  or  iii  are  of  course  omitted. 

The  variables  s,  z  and  w  are  related  as  follows  (Franklin  and  Powell 

(1980)): 

2  z-1  2  sT 

w  =  —  -  =  —  tanh  —  (4.5) 

T  z+ 1  T  2 

where  T  is  the  sampling  time,  and  the  discrete  frequecy  v  in  the  w-plane  and 

the  continuous  frequency  co  in  the  s-plane  are  related  by 
2  coT 

v  =  - tan  -  (4.6) 

T  2 

This  means,  that  for  small  sT/2  or  small  coT/2  the  following  approximations 

w  =  s  and  y  ~  lo  (4.7) 

are  valid.  Furthermore  if  G(z)  is  calculated  from  GfsJG^js)  using  the  Tustin  bilinear 
substitution  (Franklin  and  Powell  (1980)) 

2  z-1 

s  = -  (4.8) 

T  z+1 

and  G(w)  is  calculated  from  G(z)  using  the  bilinear  transformation  given  in 
equation  4.5,  then  s  =  w,  and  hence  G(w)  =  G(s)GK(s).  Therefore  evaluating 
G(s)G^(s)  for  s  =  ioo  ,  0  <  oo  <  (  gives  approximately  the  same 

Nyquist  loci  as  evaluating  G(w)  for  w  =  jv,  0  <  v  <  v0 .  The  approximation 
involved  is  valid  only  if  ooT  is  small,  which  in  practice  means  the  sampling 
time  should  be  small  relative  to  the  smallest  open  loop  time  constant.  Franklin 
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Table  4.1  Summary  of  four  different  approaches  to  the  design  of  a  discrete 
control  system  for  a  continuous  plant. 


Model 

used 

Design  and  analysis 
plane 

Quality  of  approach 

G(s) 

s-plane 

bad 

G(z) 

z-plane 

good 

G(w) 

w-plane 

good 

G(s)GK(s) 

s-plane 

fair 

and  Powell  (1980)  recommend  a  sampling  frequency  of  10  -  20  times  the 
closed  loop  bandwidth,  which  will  generally  give  a  sufficiently  small  sampling 
time  for  the  approximation  to  be  valid.  So  when  a  continuous  model  is  available 
the  fourth  approach  to  control  system  design  outlined  above  is  recommended. 

The  design  examples  of  the  next  chapter  will  use  the  third  and  fourth 
above  mentioned  approaches  to  control  system  design. 

4.3.2  Matrix  dominance  of  the  return  difference  matrix. 

The  determination  of  stability  using  the  Nyquist  exact  loci  requires  (cf. 

section  4.3.4),  that  the  return  difference  matrix  is  matrix  dominant  (cf.  chapter 
3).  As  an  initial  step  in  the  NEL  design  procedure  the  return  difference  matrix 

F(s)  =  l+G(s),  i.e.  Kj(s)  =  K^(s)  =  I,  is  tested  for  matrix  dominance  as  described 

in  chapter  3. 

If  the  return  difference  matrix  is  matrix  dominant,  the  NEL  design  can 
proceed  with  the  design  of  a  diagonal  controller  matrix  Kz(s)  as  discussed  in 
the  following  subsections,  or  a  compensator  K^(s)  can  be  designed  to  further 

reduce  interactions  or  satisfy  other  design  specifications 
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If  the  return  difference  matrix  is  not  matrix  dominant,  the  magnitude  of 
the  gains  kj(s)  in  Kz(s)  =  diag{kj(s)}  could  be  reduced,  and/or  a  compensator 
K^(s)  could  be  designed  to  make  F(s)  matrix  dominant  (Note:  In  the  limit  k;(s)-*0 
,  i  =  1 ,...., m,  F(s)  will  be  matrix  dominant  for  any  plant  G(s)). 

The  design  of  the  compensator  (s)  is  done  by  applying  either  the 
direct  Nyquist  array  technique  using  column  operations  or  pseudodiagonalization, 
the  inverse  Nyquist  array  technique  using  row  operations  or  pseudodiagonaliza¬ 
tion,  or  the  characteristic  loci  technique  using  eigenvector  alignment.  These 

techniques  are  described  in  detail  by  Kuon  (1975),  who  also  applied  them  to 
the  design  of  a  control  system  for  a  pilot  plant  evaporator  and  compared  the 

resulting  controllers.  Kuon  (1975)  found  the  three  techniques  to  give  very  similar 
results  if  the  main  design  criteria  was  to  minimize  closed  loop  interaction. 

4.3.3  Choice  of  initial  set  of  loop  gains. 

The  initial  choice  of  gains  k;(s)  in  Ka(s)  =  diag^kj(s)}  can  be  done  in 

various  ways,  eg.  arbitrarily  choose  kj(s)  =  1,  i  =  1, . m.  However,  one  practical 

approach  is  to  make  Q(s)  as  column  dominant  as  possible  by  a  diagonal 

similarity  transformation  D(s),  as  discussed  in  chapter  3.  Then  plot  the  direct 

—  1 

Nyquist  array  of  D(s)  Q(s)D(s)  with  Gershgorin  circles  superimposed  on  the 

diagonal  elements,  and  choose  a  set  of  initial  gains,  which  lie  outside  the  band 
swept  out  by  the  Gershgorin  circles.  It  should  be  noted,  that  the  Nyquist  exact 
loci  of  hj(s)  lie  inside  the  Gershgorin  band  on  q;;(s),  when  F(s)  is  dominant.  That 
is  the  magnitude  of  hj(s)  is  bounded  by: 

(Iqjj(s)l  -  Ra(s))  lh;(s)l  <  (lqfi(s)  I  +  Re(s))  (4.9) 

where,  if  F(s)  is  dominant,  Ra(s)  could  be  the  Gershgorin  radius: 

m 

R-(s)  =  H  lqu(s)l  (4.10) 

B  j=1  J 

j*i 

It  is  possible  to  define  even  narrower  bands  for  h  (s)  than  the  Gershgorin  band, 

whenever  the  cofactor  j  <s)  contains  a  row  or  column  of  zeros  Such  bands, 

called  G°-circles,  are  discussed  by  Kuon  (1975).  These  bands  of  course  have 

the  advantage  of  being  independent  of  the  elements  of  the  diagonal  controller 
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matrix  K^(s).  However,  except  for  determining  an  initial  set  of  loop  gains  bands 
are  not  used  in  the  NEL  design  procedure,  since  the  exact  loci  of  h-(s)  are 
easily  calculated  on  a  digital  computer.  In  particular,  Ostrowski  circles  are  not 
used  to  check  the  final  design  because  the  Nyquist  exact  loci  provide  a  better 
means  of  checking,  that  design  specifications  are  meet. 


4.3.4  Nyquist  exact  loci  plots. 

The  Nyquist  exact  loci  are  the  polar  plots  of  hj(s)  with  s  a  function  of 

frequency.  These  plots  are  used  in  the  NEL  procedure  to  check  if  the  system 
is  stable  and  if  the  control  is  satisfactory.  The  following  theorem  is  the  basis 

for  assessing  stability  using  the  Nyquist  exact  loci: 

Theorem  4.1:  Let  the  return  difference  matrix  F(s)  be  matrix 
dominant  and  kj(s)h j<s)  map  the  right  half  plane  Nyquist  contour  D 
into  Vh;  encircling  the  critical  point  (-1,0)  nK;  times  clockwise 

Then  the  closed  loop  system  will  be  stable  if  and  only  if 


m 


^  nkf  =  -P 


(4.1  1) 


where  pe  is  the  number  of  right  half  plane  poles  of  the  open 
loop  system. 

Proof  of  theorem  4.1  is  given  in  appendix  E. 

Theoretically  the  transfer  functions  k;(s)h.(s)  should  be  evaluated  along  the 
complete  Nyquist  contour,  but  in  practice  they  are  evaluated  only  along  a  finite 
part  of  the  imaginary  axis  (finite  frequency  range)  and  the  contour  completed  by 
extrapolation.  Also  if  kj(s)h-(s)  is  bounded  above  a  certain  frequency  the  finite 
frequency  range  results  of  Mee  (1976)  and  Leinmger  (1979b)  apply. 


4.3.5  Specify  desired  stability  margins. 

Based  on  the  Nyquist  exact  loci,  a  desired  stability  margin  is  specified 
for  each  loop.  This  margin,  as  in  classical  SISO  design,  can  be  either  a  gain 
margin  or  a  phase  margin.  The  designer  chooses  the  stability  margins  according 
to  the  design  specifications,  i.e  if  tight  control  is  desired  in  a  certain  loop  a 
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small  margin  is  specified  and  vice  versa.  The  specific  values  of  the  margins  of 

course  depend  on  the  individual  application,  but  values  of  4  -  8  db  are  typical 
for  gain  margins,  and  30  -  45  degrees  for  phase  margins. 

4.3.6  Calculate  gains  to  satisfy  stability  margins. 

Controller  constants,  which  satisfy  all  specified  stability  margins,  are 
calculated  simultaneously  using  one  of  the  algorithms  described  in  section  4.4. 

These  algorithms  will  update  the  controller  constants,  so  the  location  of  the 
Nyquist  exact  loci  approach  the  desired  locations.  This  is  done  by  linearizing  the 
underlying  optimization  problem  and  solving  a  set  of  linear  equations. 

The  updating  of  the  controller  constants  continues,  each  time  solving  a 
set  of  linear  equations,  until  the  difference  between  the  desired  locations  and 
the  actual  locations  of  the  Nyquist  exact  loci  becomes  sufficiently  small. 
Between  each  updating  of  the  controller  constants  the  designer  can  judge 

whether  to  continue  or  not  by  inspecting  the  Nyquist  exact  loci  diagrams.  The 
designer  decides  when  to  stop  the  iterative  updating  of  the  controller  constants. 

4.3.7  Is  the  control  system  satisfactory. 

The  final  step  of  the  NEL  design  procedure  is  to  check,  that  the  control 

system  is  satisfactory.  This  involves  first  of  all  checking  for  stability.  If  the 
final  return  difference  matrix  is  matrix  dominant,  then  stability  can  be  judged 
from  the  Nyquist  exact  loci,  by  counting  their  encirclements  of  the  critical  point 
(-1,0).  Since  matrix  dominance  of  F(s)  is  a  sufficient  but  not  necessary  condition 
for  stability  analysis  based  on  the  Nyquist  exact  loci,  the  system  can  be  stable 
even  though  F(s)  is  not  matrix  dominant.  Therefore  if  the  return  difference 
matrix  is  not  matrix  dominant,  stability  can  be  checked  by  a  method,  which 
does  not  require  system  dominance,  e.g.  the  characteristic  loci  technique. 

Since  the  Nyquist  exact  loci  represent  the  actual  transmittances  between 

given  input  -  output  pairs  other  control  system  performance  criteria,  such  as 
bandwidth,  rise  time,  overshoot  and  settling  time  for  unit  step  change  in  input, 

can  be  evaluated  from  the  NEL  diagrams.  These  performance  criteria  can  be 
calculated  either  from  Nyquist  plots  of  h-(s)  using  Nichols  charts,  or  from  Bode 
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plots  of  lh*,(s)!/l  1  +h-(s)l.  In  a  computer  aided  implementation  of  the  NEL  design 
technique  the  last  approach  is  more  convenient. 

The  transfer  function  hj(s)  in  general  involve  ratios  of  polynomials  in  s  of 
degree  greater  than  two.  For  such  transfer  functions  it  is  not  feasible  to  use 
analytical  expressions  for  rise  time,  overshoot  and  settling  time,  but  emperical 
correlations  have  been  published,  e.g  by  Horowitz  (1963).  The  following 

expressions  are  given  by  Horowitz: 

Rise  time  tp: 

tr  =  K  oob  /2T  0.3  <  K  <  045  (4.12) 

Overshoot  ^ 

J  =  58(  6Db  /  coA  )M  -  39  (4.13) 

Settling  time  ts< : 

ts1  =  (  ooa/4  7?  )(7(  coh!  c )M  -  3)  (4.14) 

where  40^  is  the  bandwidth  frequency,  that  is  the  frequency  at  which 

lhj(s)l/l  1  +hj(s)l  =  0.707,  and  is  the  frequency  at  which  lhj(s)l/l  1 +h;(s)l  =  0.5. 

M  is  the  magnitude  of  lhj(s)l/l  1  +hj(s)l  at  the  resonance  frequency. 

The  above  expressions  according  to  Horowitz  are  only  accurate  to  within 
15%  -  30%,  but  they  nevertheless  give  the  significant  features  of  the  time 

response  and  are  useful  in  comparing  -designs 

If  the  control  system  design  is  unsatisfactory  the  designer  can  either 

relax  the  specification  with  respect  to  stability  margins  or  other  criteria  or 

improve  the  design  of  the  compensator  K^ls)  to  meet  the  design  criteria  The 

latter  will  generally  mean  a  more  complicated  compensator. 


4.4  Simultaneous  controller  constant  evaluation. 

The  simultaneous  design  of  m  single  loop  controllers  k*(s),  i  =  1 ,m, 

each  having  n  controller  constants,  to  satisfy  desired  stability  margins  for  the 

exact  loci  of  kj(s)h;(s),  i  =  1, m,  is  mathematically  a  m  x  n  order  nonlinear 

optimization  problem.  Algorithms  to  iteratively  solve  this  problem  for  n  =  1  or 
2  are  outlined  below.  These  algorithms  use  the  Nyquist  exact  loci  and  a 
standard  Gaussian  elimination  procedure  for  adjusting  the  controller  constants  to 
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give  designer  specified  stability  margins. 

The  central  step  in  the  algorithms  is  the  solution  of  a  set  of  linear 
equations: 

A  x  =  b  (4.15) 

where  A  is  a  Jacobian  matrix  containing  the  derivatives  of  the  Nyquist  exact 
loci  w.r.t.  the  controller  gains,  x  is  a  vector  containing  the  necessary 
adjustments  of  the  controller  constants  to  satisfy  the  stability  margins,  and  b  is 
a  vector  of  distances  the  loci  have  to  be  moved  to  give  the  desired  stability 
margins,  i.e.  bj  is  the  change  in  location  of  k-,(s)h5(s).  The  elements  of  A  and  b 

are  evaluated  at  the  set  of  frequencies  s  =  JtOj  ,  i  =  1 ,m,  at  which  the 

current  gain  or  phase  margins  are  calculated.  The  distances  between  the  current 
locations  and  the  desired  locations  of  the  Nyquist  exact  loci  are  specified  by 
the  designer  via  specification  of  stability  margins.  The  actual  evaluation  of  the 
elements  of  A  and  b  depend  on  whether  n  =  1  or  n  =  2. 

For  n  =  1,  proprotional  control,  the  elements  of  A  =  a  (s)  are: 

a-j(s)  =  0  (lk;h;(s)l)/5  kj  (4.16) 

These  real  derivatives  are  of  course  evaluated  numerically  using  a  difference 
formula  For  phase  margins  the  distances  contained  in  the  vector  b  become  b; 
=  1  -  Cj,  where  c;  is  the  current  magnitude  of  kjhj(s)  at  a  phase  of  (180 

phase  margin).  For  gain  margins,  expressed  as  the  factor  the  gain  should  be 
multiplied  by  for  the  loop  to  become  unstable,  the  distances  in  the  vector  b 
become  bj  =  (1/gain  margin)  -  c;,  where  C;  here  is  the  current  magnitude  of 

k;hj(s)  at  a  location  closest  to  the  critical  point  (-1,0).  The  proportional  control¬ 
ler  constants  are  updated  using 

kp,=  kfi  +  cC  *x;  (4.17) 

where  cl  is  a  relaxation  factor  less  than  one.  The  relaxation  factor  is 
introduced  to  avoid  going  too  far,  and  possibly  oscillating  around  the  desired 
location  because  the  linearity  assumed  in  calculating  the  derivatives  a^(s),  might 
not  be  valid  for  the  full  change  in  the  controller  constants. 

For  n  =  2,  proportional  plus  integral  (PI)  or  proportional  plus  derivative 
(PD)  control,  the  elements  of  A  =  are: 

a-.(s)  =  Q  (kj(s)hj(s))/ 3  k-(s) 


(4  18) 
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and  again  these  complex  derivatives  are  evaluated  numerically  using  difference 

formula  For  phase  margins  0;  ,  i  =  1, . m,  the  elements  of  the  vector  b 

become  bj  =  (-cos  ©;  -j*sin  Q}  )  -  C;,  where  c;  is  the  location  in  the 
complex  plane  on  the  locus  kj(s)hf(s),  which  has  a  phase  of  (180  -  ©f  ).  For 

gain  margins  the  elements  of  the  vector  b  become  b j  =  (-1/gain  margin  +  j*0) 
-  Cj,  where  Cj  is  the  location  in  the  complex  plane  on  the  locus  kj(s)hj(s), 
which  is  closest  to  the  critical  point.  The  updating  of  the  controller  constants  is 
dependent  on  whether  PI  or  PD  controller  constants  are  being  calculated.  For  a 
PD-controller  of  the  form 

Kp(1+K^s)  (4.19) 

the  gains  are  updated  as  follows: 

kPi  =  kPi  +  *Re  x;  (4.20) 

and 

kDi  =  +  *(lm  xi/(Re  xi*  ))  (4.21) 

where  is  the  frequency  at  which  bj  is  calculated.  For  a  Pl-controller  of 

the  form 

Kpd+Kj/s)  (4.22) 

the  controller  constants  are  updated  as  follows: 

k-  -  k*>.  +  *Re  X.  (4.23) 

and 

kT.  =  kr.  +  *(lm  x-*  ODj  /Re  x-.)  (4.24) 

In  the  above  equations  Re  Xj  and  Im  Xj  of  course  denote  respectively  the  real 

and  imaginary  part  of  Xj. 

The  magnitudes  of  the  elements  of  the  matrix  A  provide  measures  of 

the  sensitivity  of  the  Nyquist  exact  loci  to  changes  in  any  of  the  controller 

constants.  Such  sensitivity  coefficients  can  be  a  helpful  guide  in  field  tuning  of 
a  multivariable  control  system.  These  coefficients  are  also  an  indirect  measure 
of  the  sensitivity  to  the  propagation  of  interactions.  If  the  i'th  Nyquist  exact 

locus  has  low  sensitivity  to  changes  in  the  j'th  controller  constant,  then  it  will 

likely  also  have  low  sensitivity  to,  for  example,  a  setpoint  change  in  the  j'th 

loop.  Hence  the  matrix  A  gives  a  qualitative  indication  of  which  interactions  are 

most  likely  to  upset  the  system.  The  use  of  the  matrix  A  as  a  sensitivity 


. 
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measure  warrants  further  investigation. 


4.5  Conclusions. 

The  Nyquist  exact  loci  and  the  transfer  functions  h.(s)  were  introduced. 
The  transfer  functions  h*(s)  were  shown  to  be  an  exact  multivariable  parallel  to 
the  single  variable  transfer  function  g(s).  The  Nyquist  exact  loci  were 
demonstrated  to  be  useful  in  stability  analysis  and  to  allow  control  system 

design  criteria,  such  as  rise  time,  overshoot  and  settling  time,  familiar  from 

SISO  design,  to  be  used  meaningfully  in  a  multivariable  context. 

The  Nyquist  exact  loci  (NEL)  design  procedure  was  outlined  and  its 
different  steps  discussed  in  detaii.  The  NEL  procedure  allows  the  designer  to 
specify  stability  margins  in  familiar  terms,  such  as  gain  or  phase  margin,  and 
simultaneously  design  m  single  loop  controllers.  Algorithms  for  the  simultaneous 

design  of  P,  PI  or  PD  controllers  to  give  Nyquist  exact  loci,  which  satisfy 

designer  specified  stability  margins,  were  developed  and  described 

Finally  the  use  of  the  derivatives  of  the  Nyquist  exact  loci  as  a 
sensitivity  measure  was  discussed. 


' 


5.  Control  system  design  using  NEL. 


5.1  Introduction. 

The  final  test  of  any  control  system  design  technique  is  its  use  in  a 
commercial  environment,  such  as  a  chemical  plant.  However,  before  this  can 
happen  a  design  technique’s  practicality  must  be  demonstrated  by  applying  it  to 
the  design  of  several  control  systems  for  widely  different  chemical  plants, 
comparing  the  resultant  designs  with  control  systems  designed  using  other 

techniques,  and  testing  the  designed  systems  by  simulation  and/or  application  to 
pilot  plants.  The  purpose  of  this  chapter  is  to  carry  out  this  first  part  of  the 
test  of  the  NEL  design  procedure. 

The  NEL  design  technique  will  be  used  to  design  control  systems  for 

three  chemical  processes:  a  double  effect  evaporator,  a  chemical  reactor,  and  a 
distillation  column.  These  three  processes  have  been  the  subject  of  numerous 

control  studies  in  recent  years.  This  allows  the  control  systems  designed  using 
the  NEL  procedure  to  be  compared  with  published  control  systems,  which  were 
arrived  at  using  other  design  techniques,  such  as  the  inverse  Nyquist  array 
method,  the  characteristic  loci  method,  or  the  direct  Nyquist  array  method. 

In  the  next  section  the  NEL  design  procedure  is  used  to  design  a 
discrete  control  system  for  a  pilot  plant  double  effect  evaporator.  The  design  is 
approaches  using  both  the  discrete  TFM  G(w)  and  the  continuous  TFM  G(s)G^(s), 
and  the  resulting  control  systems  are  compared  The  resulting  designs  are  also 
compared  with  the  controllers  designed  by  Kuon  (1975).  This  is  followed  by 
application  of  the  NEL  design  technique  to  the  design  of  a  control  system  for 

an  open  loop  unstable  chemical  reactor.  Finally  the  NEL  procedure  is  used  to 
design  a  control  system  for  a  distillation  column  with  significant  time  delays. 


5.2  Double  effect  evaporator. 

The  double  effect  evaporator  pilot  plant  in  the  Department  of  Chemical 
Engineering  at  the  University  of  Alberta  has  been  the  object  of  several  control 
studies,  which  were  collected  into  a  case  study  by  Fisher  and  Seborg  (1976) 
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These  studies  involved  techniques  such  as  multiloop  control,  optimal  feedback 

control,  time  optimal  control,  model  reference  adaptive  control  and  multivariable 
Smith  predictors.  Frequency  domain  control  studies  involving  the  evaporator  have 
been  reported  by  Kuor.  (1975),  who  used  the  INA,  DNA  and  CL  techniques.  Due 
to  these  past  studies  the  evaporator  is  a  good  candidate  for  applying  a  new 
design  technique,  such  as  the  NEL  design  procedure  described  in  chapter  4,  and 

comparing  the  resultant  designs  with  control  systems  obtained  using  other 
methods. 

Fisher  and  Seborg  (1976)  list  several  state  space  models  of  the 

evaporator.  Among  them  are  a  continuous  fifth  order  state  space  model  and  a 
discrete  fifth  order  state  space  model  using  a  sampling  interval  of  6A  seconds. 

The  continuous  transfer  function  matrix  (TFM!  G(s)  calculated  from  the  continuous 
state  space  model  is  shown  in  table  3.1  (Note:  all  coefficients  have  been 
truncated  to  two  significant  figures).  The  discrete  TFM  G(z)  obtained  from  the 

discrete  state  space  model  is  listed  in  table  5.1,  and  the  corresponding  bilinearly 

transformed  TFM  G(w)  is  given  in  table  5.2.  The  outputs  are  product 
concentration,  first  effect  holdup  and  second  effect  holdup,  and  the  inputs  or 

manipulated  variables  are  steam  flowrate  to  first  effect,  bottoms  flowrate  from 
first  effect  and  bottoms  flowrate  from  second  effect. 

The  poles  of  the  TFMs  G(s),  G(z)  and  G(w),  i.e.  the  roots  of  the  open 

loop  characteristic  polynomials,  are  tabled  in  table  5.3.  The  evaporator  models 

are  marginally  stable,  since  G(s)  and  G(w)  have  two  poles  located  on  the 
imaginary  axis,  and  G(z)  has  the  corresponding  two  poles  located  on  the  unit 
circle.  The  number  of  clockwise  encirclements  required  by  Nyquist  type  stability 

theorems  hence  depends  on  the  particular  Nyquist  contour  employed,  as 

discussed  in  appendix  D.  In  the  following  analysis  the  closed  right  half  plane 
Nyquist  contour  will  be  used,  this  means  two  anticlockwise  encirclements  of  the 

critical  point,  due  to  the  poles  at  the  origin  in  s-plane  or  the  w-plane. 

A  graphical  test  of  F(s)  =  l+G(s)  for  column  matrix  dominance  is  shown 

in  figure  3.1,  and  in  figure  3.2  the  corresponding  test  for  row  matrix 

dominance  is  shown  In  figure  3.1  N^  ,  N^  and  N^ ,  as  defined  by  equation 

3.6,  are  plotted  as  functions  of  frequency  for  frequencies  in  the  range  from 
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Table  5.1  Discrete  transfer  function  model  G(z)  of  pilot  plant  evaporator.  G(z)  = 


N(z)/d(z),  where  N(z)  is  a  polynomial  matrix  and  d(z)  is  the 
characteristic  polynomial. 

d(z)  = 

-0. 039+2. 49z-6  1  3z*+7.35za  -4.32z4  +  1.0zs 

n  11 <z) 

=  -0.0093+0.0  17z+0.0078za-0.029za +0.0  14z  4 

nia,<2> 

=  -0.0  1 7+0.094z-0. 1  8z*+0. 1 5z*  -0.043z4 

n,a(z) 

=  0.0 

niM  121 

=  0.0080-0.0  1 5z-0.0068z*  +  0.025z3  -0.0  1  2zs 

ntol2> 

=  -0.032  +  0. 17z-0.33z2'+0.27z3  -0.082z4 

na3<2) 

=  0.0 

n31  |Z) 

=  0.0093-0.0  1  7z-0.0079za +0.029z3  -0.0  14z4 

nM(Zl 

=  0.033-0. 1  8z+0.34z*-0.28z3  +0.085zM 

n  (z) 

3V 

=  -0.0  16+0.085z-0.16za+0.13z3 -0.041  ZH 

. 
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Table  5.2  Bilinear  transformed  discrete  transfer  function  model  G(w)  of  pilot 
plant  evaporator.  G(w)  =  N(w)/d(w),  where  N(w)  is  a  polynomial 

matrix  and  d(w)  is  the  characteristic  polynomial. 

d(w)  =  0.00005w2'  +  0.0070w3+0.24wH  +  1.0w^ 

n  1<t  (w)  =  0.00007w‘2'+0.0023w3-0.0039wV-0.0011w<s 

n,a(w)  =  -0.00005wa'-0.0026w3-0.0067w<f  +0.022w  S 

nia(w)  =  0.0 

nai(w)  =  -0.00007wa-0.0020w3 +0.0034wM +0.00098ws 
n^iw)  =  -0.00015wa-0.0050w3-0.012wV+0.041ws 
naa  (w>  ~  0.0 

n31(w)  =  -0.00008w2'-0.0023w3  +  0.0039w<<  +0.00  1  1  ws 
n  <w)  =  0.000  16wa  +  0.005  1w*+0.0  12ws -0.042ws 
n33(w)  =  -0.00008w*-0.0025w:i-0.0059wS-0.020w  s 


Table  5.3  Poles  of  evaporator  TFMs  G(s),  G(z)  and  G(w)  obtained  from  fifth 
order  state  models. 


s-plane 

z-plane 

w-plane 

0.0 

1.0 

0.0 

0.0 

1.0 

0.0 

-0.0380 

0.9603 

-0.0380 

-0.0766 

0.9215 

-0  0766 

-0.773 

0.4384 

-0.732 

0.0038  to  7.73  rad/min.  N  ^  and  are  zero  at  all  frequencies,  but  > 

1  for  frequencies  less  than  0  08  rad/min.  Hence  the  return  difference  matrix  is 
not  column  matrix  dominant  at  low  frequencies.  Similarly  in  figure  3.2  I\|J1, 
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and  as  defined  by  equation  3.6,  are  plotted  as  functions  of  frequency. 

p*  y->  f0 

N  ,a  is  less  than  one  at  all  frequencies,  but  N<3  and  are  greater  than  one 

for  frequencies  less  than  0.08  rad/min.  Therefore  the  return  difference  matrix  is 
also  not  row  matrix  dominant  at  low  frequencies.  The  return  difference  matrix 
could  be  made  matrix  dominant  by  selecting  gains  k1  and  in  diag^k1  ,k^,k3} 
less  than  one,  for  example  with  k,  =  k^  =  0.33  and  k3  =  1.0  the  return 
difference  matrix  is  matrix  dominant.  However,  such  a  scheme  would  give  loose 
control  over  the  most  important  output  variable,  i.e.  product  concentration,  which 
is  clearly  not  desirable.  The  design  of  a  compensator  to  reduce  interactions  is 
thus  warrented.  Figure  3.1  indicates,  that  the  offdiagonal  elements  of  columns  1 
and/or  2  should  be  reduced,  where  as  figure  3.2  indicates,  that  the  offdiagonal 
elements  of  row  3  should  be  reduced  to  lower  the  interaction.  This  combined 
with  a  design  objective  of  tight  control  over  product  concentration  points  to 
reducing  element  (1,2).  From  the  DNA  of  the  TFM  G(s),  which  is  shown  in 
figure  3.3,  it  appears,  that  the  open  loop  TFM  will  become  lower  triangular  by 
the  following  column  operation 

column  2  =  column  2  +  0.73*column  1  (5. 1  ) 

As  discussed  in  chapter  3  this  operation  will  tend  to  cancel  the  influence  of 
element  (2,1)  and  at  the  same  time  increase  the  direct  transmittance  in  the 

second  loop  and  slightly  decrease  the  influence  of  element  (3,2).  With  the 
above  compensation  the  return  difference  matrix  is  row  matrix  dominant,  as 
shown  in  figure  34.  The  lower  triangular  open  loop  transfer  function  matrix, 

whose  direct  Nyquist  array  is  shown  in  figure  5.1,  means  the  closed  loop 

transfer  function  matrix  will  also  be  lower  triangular,  hence  product 

concentration  is  independent  of  the  first  and  second  effect  bottoms  flows  and 
only  dependent  on  steam  flow  to  the  first  effect. 

The  simplest  compensator,  which  decouple  the  product  concentration  from 
the  bottoms  flowrates,  and  at  the  same  time  give  a  matrix  dominant  return 
difference  matrix,  seems  to  be  the  following  constant  compensator 
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Figure  5.1  DNA  of  G(s)K1  for  continuous  evaporator  TFM  G(s) 
K  i  from  equation  5.2  Frequency  range  0.0038 
rad/min. 
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1  .9 


and  compensator 
<  co  <  7.73 
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Ki  - 


1.0 

0.73 

0.0 

0.0 

1.0 

0.0 

0.0 

0.0 

1.0 

(5.2) 


If  should  be  noted,  that  G(s)K,j  is  neither  diagonally  dominant  nor  matrix 

dominant. 

The  Nyquist  exact  loci  of  G(s)K,,  G^(s),  where  GK(s)  is  a  zero  order  hold 

corresponding  to  a  sampling  time  of  64  seconds,  are  shown  in  figure  5.2. 

Careful  inspection  of  figure  5.2  shows,  that  none  of  the  loci  contribute  any 
encirclements  of  the  critical  point,  however  multiplication  of  columns  2  and  3 

by  -1  will  give  the  necessary  two  anticlockwise  encirclements  of  the  critical 
point  (-1,0),  one  each  by  qaA(s)  and  q33(s).  Hence  the  compensator  is 

TO  -0.73  0.0' 


= 


0.0  -1.0 


0.0 


0.0 


0.0 


-1.0 


(5.3) 


The  simultaneous  gain  calculation  algorithm  was  used  to  calculate 
proportional  controller  gains  to  give  specified  gain  margins.  The  results  are 
summarized  in  table  5.4.  In  all  cases  the  initial  gains  were  arbitrarily  choosen  to 
be  one,  and  after  five  iterations  the  difference  between  the  actual  and  desired 
locations  of  the  Nyquist  exact  loci  was  less  than  0.1.  Except  in  the  fourth  and 
fifth  calculation  of  table  5.4  the  return  difference  matrix  was  dominant  with  the 
final  gains.  The  two  calculations  in  which  F(s)  was  not  dominant,  are  the  ones 
with  the  largest  relative  difference  between  the  gain  margin  in  loop  1,  the 
concentration  loop,  and  those  in  loops  2  and  3,  the  level  loops  In  the  fourth 
calculation  it  is  evident  from  the  row  matrix  dominance  plot,  which  is  shown  in 
figure  5.3,  that  F(s)  can  be  made  matrix  dominant  by  a  slight  increase  in 
The  results  of  calculations  number  one  and  two,  due  to  the  large  ratio  ka/k-j, 
are  not  in  agreement  with  the  objective  of  tight  control  of  product 
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Figure  5.2  NEL  of  continuous  evaporator  with  TFM  G(s),  zero  order  hold  with 
sampling  time  64  seconds,  compensator  K,  from  5.3  and  =  I. 
Frequency  range  0  00,38  <  co  <  7.73  rad/min. 
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Fiqure  5  3  Row  matrix  dominance  test  of  continuous  evaporator  with 

compensator  K  ^  from  equation  5.3,  zero  order  hold  with  64 
second  sampling  time,  and  proportional  gain  from  calculation  no.  4 
of  table  5.4. 
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Table  5.4  Summary  of  proportional  controller  constant  calculation  for  continuous 
evaporator  model  G(s)Gk(s)K< . 


Calc.  no. 

gain 

margins 

proportional 

constants 

F(s)  dominant  ? 

1 

2 

2 

2 

10.1 

14.1 

40.5 

yes 

2 

2 

3 

3 

9.9 

5.9 

23.5 

yes 

3 

2 

3 

4 

9.9 

5.9 

16.5 

yes 

4 

2 

4 

4 

10.6 

3.4 

16.5 

no 

5 

2 

5 

5 

1  1.3 

2.3 

12.6 

no 

6 

3 

4 

4 

5.6 

4.1 

15.6 

yes 

7 

3 

5 

5 

5.8 

2.9 

1  1.8 

yes 

8 

2.5  4 

4 

7.2 

3.6 

15.6 

yes 

concentration,  and  loose  control  of  the  two  levels.  This  leaves  the  results  of 

calculations  number  six  to  eight,  in  figure  5.4  to  5.6  lhj(s)l/l  1  +hj(s)l  for  i  =  1,2,3 

are  plotted  as  functions  of  frequency.  The  proportional  gains,  which  best  fulfills 
the  objective  of  tight  control  on  product  concentration  and  loose  level  control 
are  those  of  calculation  number  eight.  The  simultaneous  gain  calculation  took 
four  steps  to  arrive  at  the  final  gain  starting  with  k1  =  ka  =  k3  =  1.  The 

Nyquist  exact  loci  after  each  iteration  are  shown  in  figures  5.7  to  5.10.  It  is 
evident  from  these  figures,  that  the  shape  of  the  Nyquist  exact  loci  do  not 

change  drastically  from  iteration  to  iteration,  and  figure  5.10  also  shows  the 
desired  gain  margins  have  been  realized.  Rise  time,  overshoot  and  settling  time 

as  calculated  using  equations  4.12  to  4.14  with  the  gains  of  calculation  number 

eight  are  shown  in  table  5.5  Loop  1,  the  product  concentration  loop,  responds 
purest  to  a  unit  setpoint  change  However,  since  all  variables  are  normalized 

deviation  variables,  a  unit  step  in  product  concentration  actually  represents  a 
doubling  of  the  product  concentration,  which  is  physically  more  difficult  to 
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Figure  5.4  Plots  of  lk-t h  (s)l/l  1  +k-,hj(s)l  versus  frequency  for  i=  1,2,3  with 

proportional  gains  from  calculation  no.  6  of  table  5.4. 


90 


Figure  5.5  Plots  of  lkjh-((s)l/N +kjh;(s)l  versus  frequency  for  i=  1,2,3  with 
proportional  gains  from  calculation  no.  7  of  table  5.4. 
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Figure  5.6  Plots  of  lkjhj(s)l/l  1  +k;h5(s)l  versus  frequency  for  i=  1,2,3  with 

proportional  gains  from  calculation  no.  8  of  table  5  4. 
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Figure  5.7  Nyquist  exact  loci  for  continuous  evaporator  model  after  one 
iteration  of  gain  calculation  no.  8  of  table  5.4 


. 
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Figure  5.8  Nyquist  exact  loci  for  continuous  evaporator  model  after  two 
iteration  in  gain  calculation  no.  8  of  table  5.4. 


■ 
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Figure  5.9  Nyquist  exact  loci  for  continuous  evaporator  model  after 
iterations  in  gain  calculation  no.  8  of  table  5.4. 


three 
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Fiaure  5.10  Nyquist  exact  loci  for  continuous  evaporator  model  after  last 
(fourth)  iteration  in  gain  calculation  no.  8  of  table  5.4. 


. 
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Table  5.5  Rise  time,  overshoot  and  settling  time  with  proportional  gain  from 
calculation  number  eight  of  table  5.4. 


loop  1 

loop  2 

loop  3 

rise  time,  min. 

4.7 

3.4 

1.9 

overshoot, % 

1  1 

5 

3 

settling  time,  min. 

12.5 

8.7 

4.5 

accomplish  than  a  doubling  of  the  hold-ups  of  each  effect  The  recommended 
final  proportional  compensator/controller  hence  is 


7.2 

-2.6 

0.0 

K-i  Ka  - 

0.0 

-3.6 

0.0 

_  0.0 

0.0 

-15.6  . 

This  compensator/controller  can  be  compared  with  the  proportional  controllers 
FD0320  and  FD0330  designed  by  Kuon  (1975): 


CO 

ib 

i _ 

-3.1 

0.0 

FD0320  = 

-2.2 

-3.7 

0.0 

(5.5) 

-9.8 

-4.9 

-9.8. 

and 

8.3 

-4.7 

o.o" 

FD0330  = 

-3.3 

-5.5 

0.0 

(5.6) 

14.8 

-7.3 

- 1 4.8_ 

which 

were  designed  for 

gain 

margins 

of  respectively  5  and  3.3 

in  all  three 

loops.  The  additional  non-zero  elements  below  the  diagonal  in  the 
compensator/controllers  FD0320  and  FD0330  will  reduce  level  oscillations,  but 
they  will  not  influence  the  product  concentration  response  to  setpoint  changes. 
Kuon  used  the  discrete  evaporator  model  G(w)  in  designing  FD0320  and 
FD0330,  and  the  controllers  were  found  to  perform  satisfactorily  in 
experimental  test.  Compared  with  Kuon’s  compensator/controllers  K<,K^  seems  to 
have  lower  gains  in  loops  1  and  2  and  higher  in  loop  3,  when  the  gain 
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margins  used  in  calculating  are  considered.  This  could  be  due  to  the  use 

of  the  continuous  model  G(s)G^(s)  in  calculating  K i .  In  order  to  investigate  this 
possibility  the  calculations  are  repeated  using  the  discrete  model  G(w). 

The  DNA  of  the  bilinear  transformed  discrete  evaporator  model  of  table 
5.2  is  shown  in  figure  5.11  for  discrete  frequencies  from  0.01  to  2.80 
rad/min.,  corresponding  to  continuous  frequencies  from  0.01  to  7.73  rad/min.  A 
compensator  identical  to  the  one  in  equation  5.2  makes  the  return  difference 
matrix  dominant  with  =  I.  Proportional  controller  gains  were  calculated  for 
the  same  conditions  as  were  used  with  the  continuous  model.  The  results  of 
the  eight  calculations  with  the  discrete  model  are  summarized  in  table  5.6.  With 
respect  to  stability  the  results  are  identical,  the  gains  of  loop  1  and  2  are 
very  similar,  but  the  gain  in  loop  3  is  consistently  lower  when  the  discrete 
model  is  used.  This  points  to  an  error  in  either  the  discrete  or  continuous 
model  state  space  model  given  by  Fisher  and  Seborg  (1976).  The  results  in 
table  5.6  of  course  agree  with  the  compensator/controllers  FD0320  and 
FD0330  (shown  in  equations  5.5  and  5.6)  designed  by  Kuon  (1975)  using  a 

discrete  evaporator  model. 

Comparison  of  tables  5.4  and  5.6  leads  to  the  conclusion,  that  the 

continuous  TFM  G(s)G^(s)  can  safely  be  used  in  place  of  the  bilinear  transformed 
discrete  TFM  G(w)  when  designing  a  discrete  control  system  for  a  continuous 
plant,  even  though  the  approximations  of  equation  4.7  are  invalid.  The  use  of 
the  continuous  TFM  G(s)G^(s)  has  the  advantage,  although  not  explored  in  the 

above  example,  that  the  effect  of  different  sample  times  can  easily  be 

investigated  without  time  consumming  transformations. 

The  simulated  closed  loop  responses  of  the  evaporator  with  the 
compensator/controller  of  equation  5.4  to  step  changes  in  setpoint  of  +10%  in 
product  concentration,  y1?  and  +20%  in  first  effect  hold-up,  ya,  are  plotted  in 
figures  5.12  and  5.13  respectively  The  discrete  evaporator  fifth  order  state 
space  model  was  used  in  these  simulations.  It  can  be  seen,  that  the  objective 
of  minimizing  effects  of  level  changes  on  product  concentration  has  been 
achieved.  These  results  compare  favourably  with  those  from  other  multivariable 
proportional  controllers  designed  for  the  evaporator,  Fisher  and  Seborg  (1976), 
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FjqUre  5.11  DNA  of  bilinear  transformed  discrete  evaporator  model  G(w). 
Discrete  frequency  range  0.01  <  v  2.80  rad/min. 
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Figure  5.12  Simulated  response  of  evaporator  to  a  +10%  change  in  product 
concentration.  The  fifth  order  discrete  state  space  model  and  the 
compensator/controller  of  5.3  were  used 
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Figure  5  13  Simulated  response  of  evaporator  to  a  +20%  change  in  first  effect 
hold-up  The  fifth  order  discrete  state  space  model  and  the 
compensator/controller  of  equation  5.3  were  used 
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Table  5.6  Summary  of  proportional  controller  constant  calculation  for  discrete 
evaporator  model  G(w)K^ . 


Calc.  no. 

gain 

margins 

proportional 

constants 

F(s)  dominant  ? 

1 

2 

2 

2 

1  1.2 

12.2 

25.1 

yes 

2 

2 

3 

3 

9.9 

7.8 

16.6 

yes 

3 

2 

3 

4 

9.9 

7.8 

12.4 

yes 

4 

2 

4 

4 

1  1.1 

3.6 

12.4 

no 

5 

2 

5 

5 

1  1.8 

2.5 

9.9 

no 

6 

3 

4 

4 

5.9 

4.5 

12.4 

yes 

7 

3 

5 

5 

6.1 

3.2 

9.9 

yes 

8 

2. 

5  4 

4 

7.5 

4.1 

12.4 

yes 

including  those  designed  using  other  frequency  domain  techniques,  Kuon  (1975). 
The  rise  time,  overshoot  and  settling  time  of  the  simulation  responses  are 
qualitatively  in  agreement  with  the  values  given  in  table  5.5. 
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5.3  Chemical  reactor. 

A  model  of  an  open  loop  unstable  chemical  reactor  has  been  used  in 

testing  control  system  design  techniques  by  Belletrutti  (1972),  MacFarlane  and 
Kouvaritakis  (1977),  and  Hung  and  Anderson  (1979).  The  different  authors  give 
slightly  different  continuous  TFM  models  of  the  reactor.  However,  all  the  models 
have  two  inputs,  two  outputs  and  two  right  half  plane  poles,  which  makes  the 
reactor  an  open  loop  unstable  system.  In  this  design  exercise  the  model  used 
by  Hung  and  Anderson  (1979)  and  listed  in  table  5.7  will  be  employed.  This 

model  has  the  right  half  plane  poles  at  0.0974  and  1.77  respectively. 

The  return  difference  matrix  F(s)  =  l+G(s)  is  matrix  dominant  as  shown  in 

figure  5.14,  where  N1z(s)  is  plotted  as  a  function  of  frequency.  N12t(s)  is  less 

than  one  at  all  frequencies,  but  the  diagonal  elements  of  the  DNA  of  G(s)  give 

only  one  anticlockwise  encirclement  of  the  critical  point,  as  evident  from  figure 
5.15.  Table  5.7  and  figure  5.15  also  show,  that  the  elements  of  the  first 

column  of  G(s)  are  two  orders  of  magnitude  smaller  than  the  elements  of  the 

second  column.  Multiplying  the  first  column  of  G(s)  by  -100  gives  the  DNA 

displayed  in  figure  5.16.  Now  the  diagonal  elements  give  a  total  of  two 
anticlockwise  encirclements  of  the  critical  point,  but  the  return  difference  matrix 
is  no  longer  matrix  dominant.  Figure  5.16  suggests  element  (2,1)  can  be  reduced 

by  the  column  operation 

column  1  -  column  1  -  1.6*column  2  (5.7) 

The  return  difference  matrix  F(s)  =  l+G(s)K  with 


-100.0 

0.0 

-1.6 

1.0 

(5.8) 


is  matrix  dominant,  see  figure  5.17,  and  the  diagonal  elements  of  the  DNA 
G(s)Ki,  gives  two  anticlockwise  encirclements  of  the  critical  point.  Hence,  the 
system  G(s)K1  is  closed  loop  asymptotically  stable  according  to  theorem  3.2. 
G(s)K1  is  matrix  dominant,  but  not  diagonally  dominant,  cf.  figure  5.18,  but  it 
can  be  made  diagonally  dominant  by  transfer  of  dominance  as  described  in 
section  3.5.  The  DNA  of  G(s)K1  after  a  diagonal  similarity  transformation 
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Figure  5.14  Matrix  dominance  test  of  F(s)  =  l+G(s)  with  the  chemical  reactor 
TFM  Gis)  from  table  5.7. 
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Fiaure  5.15  DNA  of  chemical  reactor  TFM  G(s)  from  table  5.7.  Frequency  range 
0.002  <  (jo  <  30.5  rad/ 1000s. 


105 


30.0-1 


Figure  5  16  DNA  of  chemical  reactor  TFM  G(s)  from  table  5.7  after 
multiplication  of  column  one  by  -100.  Frequency  range  0.002  < 
oo  <  30.5  rad/ 1000s. 
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Figure  5 


17  Matrix  dominance  test  of  F(s)  -  l+G(s)K1  with  the  chemical  reactor 
TFM  G(s)  from  table  5.7  and  the  compensator  K1  from  equation 

5.8. 
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53.  On 


Figure  5. 1 8  DNA  of  compensated  chemical  reactor  TFM  model  with  Gershgorin 
circles,  based  on  column  sums,  superimposed  on  the  diagonal 
elements.  Frequency  range  0.002  t o  <  30.5  rad/ 1000s. 
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Table  5.7  TFM  model  of  chemical  reactor  as  given  by  Hung  and  Anderson 
(1979).  The  model  G(s)  has  the  form  G(s)  =  N(s)/d(s),  where  d(s)  is 
the  characteristic  polynomial,  and  N(s)  is  a  polynomial  matrix. 


d(s)  =  7.55-79.43s+18.47sa  +  1  1.87s3  +  1.0sV 


n  u  (s)  =  -2.63-1  1.1  ls-2.1  Is* 


n12i(s)  =  233. 5+29. 2s 
nA1  (s)  =  1.5  15  +  0. 943s 


n^(s)  =  -96. 5-58.25  s+43. 83s*  +  5.68s"3 


calculated  with  column  weighting  factors  of  3  and  1  is  presented  in  figure 

5.19.  From  figure  5.19  it  can  be  seen,  that  the  system  will  be  stable  with 

proportional  gains  in  the  following  ranges 

0.06  4  k1  <=o  (5.9) 

0.58  <  k*  <*°  (5.10) 

The  Nyquist  exact  loci  with  k.j  =  0.15  and  k^  =  0.15  are  shown  in  figure 

5.20.  These  gains  give  gain  margins  of  23  and  3  respectively.  As  can  be  seen 

the  corresponding  phase  margins  are  60  and  90  respectively.  The  gain  is 

outside  the  range  given  in  inequality  5.10,  but  the  return  difference  matrix  is 
dominant,  as  shown  in  figure  5.21,  and  the  Nyquist  exact  loci  of  figure  5.20 
give  a  total  of  two  counterclockwise  encirclements  of  the  critical  point,  so  the 
system  is  asymptotically  stable  according  to  theorem  4.1.  Through  examination  of 
plots  of  lk*hj(s)l/l  1  +kjhj(s)l  for  different  K,^  =  diag{k,k$  it  was  found,  that 

improved  control  was  obtained  with  higher  gains.  It  was  also  observed,  that  the 
bandwidth  was  proportional  to  the  gain.  This  is  in  agreement  with  the 
controllers  designed  by  Belletrutti  (1972),  by  MacFarlane  and  Kouvaritakis  (1977), 

and  by  Hung  and  Anderson  (1979).  The  limitations  on  the  proportional  gains  will 
in  practice  be  dependent  on  the  error  in  the  linear  model  of  the  reactor.  Plots 

of  lk-,hj(s)l/l  1  +kjh-,(s)l  for  k1  =  kx  =  1.0  are  shown  in  figure  5.22,  from  which 

rise  time,  overshoot  and  settling  time  can  be  calculated  using  equations  4.12  to 
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Figure  5.19  DNA  of  similarity  transformed  compensated  chemical  reactor  TFM 
model  with  Gershgorin  circles,  based  on  column  sums,  superimposed 
on  the  diagonal  elements.  Frequency  range  0.002  <  <  30.5 

rad/ 1000s. 
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Figure  5.20  Nyquist  exact  loci  of  compensated  chemical  reactor  model  with  k, 
=  0.15  and  =  0.15.  Frequency  range  0.002  <  co  <  30.5 

rad/ 1  000s. 
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Figure  5.21  Matrix  dominance  test  of  F(s)  =  l+GislK,  with  the  chemical 
reactor  TFM  model  from  table  5.7,  the  compensator  K1  from 
equation  5.8  and  =  diag  £0. 1  5,0. 1  5} . 


. 


4.14.  The  first  loop  has  a  rise  time  of  2  min.,  65%  overshoot,  and  a  settling 
time  of  20  min.,  while  the  rise  time  of  the  second  loop  is  59  min.  The  high 
overshoot  could  be  reduced  by  introducing  integral  action  into  the 
compensator/controller.  The  final  proportional  compensator/controller  is 


K,K 


X. 


-100.0 

0.0 

-1.6 

1.0 

(5.1  1) 


It  should  be  noted,  that  the  low  gain  limits  of  the  gain  gain  space  with  the 
constant  compensator  in  equation  5. 1  1  are  somewhat  lower,  than  those  reported 
by  Belletrutti  (1972)  using  the  characteristic  locus  desing  technique  and  a  high 
integrity  diagonal  controller  matrix. 
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Figure  5.22  Plots  of  lk;hj(s)l/l  1  +k|h-,(s)l  for  compensated  chemical  reactor 
with  k,  =  k^  =  1.0.  Frequency  range  0.002  co  < 

rad/ 1  000s. 


model 

30.5 


■ 


5.4  Distillation  column. 

The  binary  pilot  scale  distillation  column  in  the  Department  of  Chemical 
Engineering  at  the  University  of  Alberta  has  also  been  the  subject  of  numerous 
control  application  studies,  e.g.  Berry  (1973)  and  Bilec  (1979). 

Berry  (1973)  obtained  a  simple  transfer  function  model  of  the  distillation 
column  by  pulse  testing.  This  transfer  function  model  G(s)  is  given  in  table  5.8. 
The  model  outputs  are  top  and  bottom  composition,  and  the  inputs  are  reflux 
flowrate  and  flowrate  of  steam  to  reboiler.  As  can  be  seen  from  figure  5.23 
the  return  difference  matrix  F(s)  =  l+G(s)  is  matrix  dominant.  It  is  not  diagonally 
dominant  at  low  frequencies,  but  G(s)  is  both  matrix  and  column  diagonally 
dominant.  However,  this  example  illustrates  how  the  NEL  design  procedure  can 
be  applied  to  systems  with  time  delays. 

The  transfer  functions  h j(s)  for  systems  with  time  delays  become  ratios 

of  polynomials  in  s  with  non-constant  coefficients.  For  example,  i\(s)  is  for 

the  distillation  column  model  in  table  5.8: 

-19.4e“3s  -ki  *1  8.9*6.6e”10*  (16.7s+1) 

h  (s)  =  -  -  - - 7—  (5.12) 

14.4s+1  (2  1.0s+  1)(  1  0.9s+  1 )( 1 6.7s+  1  +k1  *  1 2.8e”1s  ) 

The  stability  of  transfer  functions  of  this  form  has  been  studied  by  Lee  (1976), 
who  used  basic  direct  Nyquist  stability  analysis.  Hence,  the  stability  of  system 
with  time  delays  can  be  analyzed  using  the  Nyquist  exact  loci  (Note:  the  time 
delays  are  easily  handled  by  the  computer  aided  design  package,  since  the  TFM 
manipulations  of  the  NEL  design  procedure  are  implemented,  not  algebraically,  but 
numerically). 

Since  G(s)  is  column  diagonally  dominant  initial  gain  estimates  for  the 
simultaneous  calculation  of  the  two  proportional  gains  are  best  selected  using 
Gershgorin  circles  on  the  diagonal  elements  g11  (s)  and  g^(s)  of  the  DNA  of 
G(s).  The  DNA  of  G(s)K1  with  Gershgorin  circles  superimposed  on  the  diagonal 
elements  is  shown  in  figure  5.24  with  the  precompensator 
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Figure  5.23  Matrix  dominance  test  of  F(s)  -  l+G(s)  with  the  distillation  column 
TFM  G(s)  from  table  5.8. 


■ 


Table  5.8  TFM  model  of  binary  pilot  scale  distillation  column  as  qiven  by  Berry 
(1973). 

g  (s)  =  12.8e~1s/(16.7s+1) 
g1z(s)  =  -  1 8.9e""Ss/(2 1.0s+ 1) 
gz1  (s)  =  6.6e'7i/(10.9s+1) 
gzx(s)  =  -19.4e“3s/(14.4s+1) 


1.0 

0.0 

0.0 

-1.0 

(5.13) 


Since  the  distillation  column  TFM  G(s)  has  no  poles  or  zeros  in  the  closed  right 

half  plane,  the  diagonal  elements  of  G(s)K<,  should  give  zero  encirclements  of 
the  critical  point  (-1,0).  This  occurs  if  k  1  and  of  =  diag{k1  ,k^l  lie  in 
the  following  ranges 

0.0  ^  k1  <  0.89  (5.14) 

0.0  4  kx  <  0.25  (5.15) 

The  upper  limits  were  determined  from  the  band  swept  out  by  the  Gershgorin 
circles  in  figure  5.24.  Proportional  controller  constants  were  then  calculated  for 

several  gain  margin  specifications  and  initial  gains  of  0.45  and  0  12  respectively. 
The  results  of  these  calculations  are  summarized  in  table  5.9  Only  calculations 
number  2  and  3  of  table  5.9  do  not  give  large  dual  resonance  peaks  in  the 
plots  of  lkjhj(s)l/l  1  +kjhj(s)i  versus  frequency  as  seen  in  figures  5.25  and  5.26. 
The  proportional  gains  of  these  two  calculations  are  comparable  to  those  used 
by  Wood  and  Berry  (1973)  in  experimental  test  on  the  distillation  column  It  is 
not  possible  to  use  equations  4.12  to  4.14  due  to  the  nature  of  the  plots  in 
figures  5.25  and  5  26.  A  large  amount  of  interaction  exist  in  the  system,  as 
the  size  of  the  coefficients  of  the  sensitivity  matrix  A  from  calculation  number 


2  shows: 


10.0-1 


Figure  5.24  DNA  of  compensated  distillation  column  TFM  model  with  Gershgorin 
circles,  based  on  column  sums,  superimposed  on  the  diagonal 
elements.  Frequency  range  0.005  1.0  rad/s. 


' 
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Figure  5.25  Plots  of  lkjh|(s)l/l  1+kjhj(s)l  versus  frequency  for  compensated 
distillation  column  model  with  proportional  gain  from  calculation 
number  2  of  table  5.9 
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Figure  5.26  Plots  of  lk-,hj(s)l/l  1  +k;h;(s)l  versus  frequency  for  compensated 
distillation  column  model  with  proportional  gain  from  calculation 
number  3  of  table  5.9 


■ 
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Table  5.9  Summary  of  proportional  controller  constant  calculation  for 
compensated  distillation  column  model  G(s)K., . 


Calc.  no. 

gain  margins 

proportional 

constants 

F(s)  dominant  ? 

1 

2 

2 

0.89 

0.12 

yes 

2 

3 

4 

0.56 

0.085 

yes 

3 

4 

4 

0.43 

0.068 

yes 

4 

3 

2 

0.53 

0.13 

yes 

5 

4 

2 

0.40 

0.14 

yes 

A 


0.72  0.25' 

0.20  4.51 


(5.16) 


The  proportional  controllers  will  not  give  tight  composition  control.  If  tight 

control  of  both  compositions  are  required  a  dynamic  compensator,  which 
essentially  decouples  the  system,  must  be  designed. 

The  proportional  gains  of  all  calculations  fall  within  the  ranges  of 

inequalities  5  14  and  5.15,  and  it  is  hence  not  necessary  to  test  for  stability 
using  the  Nyquist  exact  loci.  However,  the  Nyquist  exact  loci  with  the  gains 
from  calculation  number  2  of  table  5.9  are  shown  in  figure  5.27.  It  can  be 

seen,  that  the  Nyquist  exact  loci  give  no  encirclements  of  the  critical  point,  so 
the  system  is  asymptotically  stable  according  to  theorem  4.1. 

An  attempt  to  use  the  algorithm  outlined  in  section  4.4  to  calculate  the 
constants  of  two  proportional  plus  integral  controllers  failed.  This  is  probably 

happening  because  the  crossover  frequency  is  not  used  in  estimating  the  integral 
gains.  The  algorithms  of  section  4.4  for  two  constant  controller  calculation 
should  be  modified  to  estimate  integral  and  derivative  gains  based  on  crossover 
frequency,  and  then  calculate  proportional  gains  to  give  desired  stability  margins. 
Such  a  modification  would  eliminate  the  use  of  complex  derivatives  and  the 
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Figure  5.27  Nyquist  exact  loci  for  compensated  distillation  column  TFM  model 
with  proportional  gain  from  calculation  number  2  of  table  5.9. 
Frequency  range  0.005  <  to  ^  1.0  rad/min. 
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mappings  of  equations  4.20  and  4.22. 


5.5  Conclusions. 

The  Nyauist  exact  loci  (NEL)  design  technique  was  applied  to  design  of 

constant  compensator/controllers  for  a  double  effect  evaporator  model,  an  open 
loop  unstable  chemical  reactor  model,  and  a  binary  distillation  column  model  with 
time  delays.  The  NEL  procedure  was  found  to  be  easy  to  use,  and  the  resultant 
constant  compensator/controllers  were  comparable  to  compensator/controllers 
designed  using  other  techniques.  The  Nyquist  exact  loci  were  found  to  be 

informative  in  comparing  design  alternatives,  and  the  absense  of 

Gershgorin/Ostrowski  circles  made  the  displays  more  comprehendable. 

The  double  effect  evaporator  design  example  showed,  that  a  continuous 
transfer  function  model  with  a  zero  order  hold  is  a  good  approximation  to  a 
bilinear  transformed  discrete  transfer  function  model.  The  constant 

compensator/controller  designed  using  NEL  was  simpler,  than  those  designed  by 
Kuon  (1975),  due  to  the  less  restrictive  matrix  dominance  condition  for  stability. 
Simulations  showed  the  compensator/controller  designed  using  NEL  performed  as 
well  as  those  designed  by  Kuon  (1975). 

The  open  loop  unstable  chemical  reactor  example  showed,  that  a  simple 
constant  compensator  gave  a  larger  stable  gain  space,  than  found  by  the 
characteristic  loci  design  procedure.  This  example  also  showed  the  estimate  of 
the  stable  gain  space  obtained  from  the  DNA  with  Gershgorm  circles  after 
transfer  of  dominance  is  conservative,  and  it  can  be  enlarged  using  the  NEL 

design  technique. 

The  binary  distillation  column  example  showed,  that  the  NEL  design 
procedure  can  be  applied  to  systems  containing  time  delays  with  no  difficulty. 
The  proportional  gain  obtained  using  the  NEL  technique  were  comparable  to 
those  used  in  experimental  test  on  the  distillation  column  by  Wood  and  Berry 
(1973).  This  example  also  showed,  that  the  algorithms  for  calculating  proportional 

plus  integral  and  proportional  plus  derivative  controllers  described  in  section  4.4 

need  to  be  modified  to  consider  the  crossover  frequency  in  the  estimation  of 


integral  and  derivative  constants. 
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6.  Conclusions  and  recommendations. 


6.1  Contributions  of  this  work. 

The  main  contributions  of  this  work  are: 

i-  Analysis  of  interaction  in  multivariable  control  systems  with  the 

introduction  of  a  classification  of  transmittances  into:  direct,  parallel, 
interaction  and  disturbance  transmittances. 

ii.  Critical  review  of  published  measures  of  'interaction  with  a 

clarification  of  the  difference  between  interaction  and  parallel 
transmittance. 

iii.  Definition  of  matrix  dominance  as  a  less  restrictive  condition  for 

the  application  of  the  multivariable  Nyquist  array  stability  theorems. 

iv.  Development  of  a  graphical  test  for  matrix  dominance,  which  also 

gives  guidance  in  the  design  of  a  compensator  for  a  system, 

which  is  not  matrix  dominant. 

v.  Development  of  a  systematic  and  flexible  procedure  for  transfer 

of  dominance  in  transfer  function  matrices  of  any  size. 

vi.  Development  of  the  Nyquist  exact  loci  (NEL)  design  techniques  as  a 

direct  multivariable  parallel  to  classical  SISO  control  system  design. 

vii.  Application  of  the  NEL  design  technique  to  three  different  systems, 

and  comparison  of  the  resultant  multivariable  controller  with  control 
systems  designed  using  other  methods. 

Although  not  a  major  contribution  this  work  has  also  involved  the  specification 

and  partly  implementation  of  a  suite  of  interactive  computer  programs  for  using 

the  NEL,  DNA  and  INA  design  procedures.  The  specification  of  the  package  is 
contained  in  Jensen  (1980).  The  programs  were  used  in  this  work,  and 

extensively  tested  during  a  multivariable  control  system  design  exercise  in  the 
course  Ch.E  644.  Students  with  no  previous  exposure  to  the  University  of 
Alberta  computing  f acidities  used  the  programs  without  difficulty  after  only  a 
short  introduction  to  the  package. 


124 


. 


125 


6.2  Further  research. 

Some  recommendations  for  future  research  in  this  area  include: 

The  following  three  recommendations  reflects  the  trend  in  multivariable 
frequency  domain  design  to  design  compensators  to  reduce  interaction  and/or 
satisty  an  abstract  condition,  like  diagonal  dominance,  and  only  indirectly  consider 
time  domain  specifications  such  as  rise  time,  overshoot  and  settling  time 

i.  Development  of  algorithms  for  simultaneous  calculation  of  controller 

constants  in  several  single  loop  controllers  subject  to  specified 
stability  margins  and  to  the  condition,  that  the  return  difference 
matrix  is  matrix  dominant.  A  linear  programming  approach  combining 
the  transfer  of  dominance  algorithm  of  section  3.5  and  the 

simultaneous  gain  calculation  of  section  4.4  should  definitely  be 
explored. 

ii.  Further  work  should  be  done  in  the  parallelism  between  the 

transfer  function  h.(s)  and  the  SISO  transfer  function  g(s),  for 
example  the  use  of  derivatives  of  the  Nyquist  exact  loci  as  an 

indicator  of  interaction  propagation  should  be  analyzed  further.  Any 
relationship  between  these  derivatives  and  the  closed  loop 

sensitivities  should  be  investigated. 

iii.  Development  of  an  algorithm  for  simultaneous  design  of  several 

lead  -  lag  elements  to  give  desired  shape  of  the  hj  loci. 

Other  recommendations  in  the  general  area  of  frequency  domain  control  system 
design  include: 

i.  Development  of  a  design  procedure  for  the  design  of  multivariable 

feedforward  controllers  for  nonlinear  plants  with  large  parameter 
variations.  The  approach  taken  by  Horowitz  (1979,1980)  to  the 
design  of  feedback  controllers  for  such  systems  should  be 
considered  for  extension. 

ii.  Analysis  of  connections  between  the  transfer  functions  h  (s)  and 


. 
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stability  criteria  for  nonlinear  systems  should  be  done,  cf.  Blight 
and  McClamroch  (1977). 

iii.  Development  of  a  graphical  representation  which  will  provide 

sufficient  conditions  for  the  Popov  criterion  to  hold  in  the 
multiloop  case  for  matrix  dominant  systems.  This  is  an  extension  to 
work  recently  reported  by  Mees  and  Atherton  (1980). 


7.  Nomenclature. 


7.1  Technical 

ADGA 

CL 

DNA 

INA 

NEL 

RDGA 

RGA 

SISO 

TFM 


abbreviations. 

Average  Dynamic  Gain  Array. 
Characteristic  Loci. 

Direct  /Vyquist  Array. 

/nverse  /Vyquist  Array. 
/Vyquist  Dxact  Loci. 

Relative  Dynamic  Gain  Array. 
Relative  Gain  Array. 

Single- input,  single-output. 
Transfer  function  Matrix. 


7.2  Nomenclature 
e 

g 

Gf> 

hi 

H 


K 

K 

m 

Q 

r 

u 

y 

y 


for  chapter  1. 

vector  of  error  signals. 

transfer  function  of  single-input,  single-output  system, 
plant  load  transfer  function  matrix, 
plant  transfer  function  matrix. 

transfer  function  between  the  i'th  input  and  the  i'th  output 

when  the  i’th  loop  is  open  and  all  other  loops  are  closed, 
feedback  transfer  function  matrix, 

indentity  matrix. 

compensator  transfer  function  matrix, 
diagonal  controller  transfer  function  matrix, 

vector  of  manipulated  varibles. 

compensated  plant  transfer  function  matrix,  Q  =  Gr,K1. 

vector  of  reference  input  variables. 

vector  of  compensated  plant  input  variables. 

vector  of  plant  output  variables. 

vector  of  measured  plant  output  variables. 

vector  of  plant  load  or  disturbance  variables. 
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7.3  Nomenclature 

A 

b; 

B 


C 

Ci 

D(t) 

F  =  Cf5il 

Gl  =  tgjil 

g 


(IE); 

H={h,p 
K,  =  {kjjl 
Kj^diagtkjl 


M!s)=!p;j(si3 
M(t>=  fyijjlt)} 
Mi-l+BjlC.Bj) 


for  chapter  2. 

state  matrix  of  state  space  model. 

the  i'th  column  of  the  state  space  model  input  matrix  B 

input  matrix  of  state  space  model. 

input  matrix  of  state  space  model  without  the  i'th  column, 

the  i'th  row  of  state  space  model  output  matrix  C. 

cofactor  of  element  (i,j)  of  the  return  difference  matrix, 
output  matrix  of  state  space  model. 

output  matrix  of  state  space  model  without  the  i'th  row. 

matrix  of  average  dynamic  gains  calculated  by  integrating  step 

responses  over  a  time  interval  t  =  t^  -  ti. 

return  difference  transfer  function  matrix. 

plant  load  transfer  function  matrix. 

plant  transfer  function  matrix. 

transfer  function  between  the  i'th  input  and  the  i'th  output 

when  the  i'th  loop  is  open  and  all  other  loops  are  closed, 

interaction  coefficient  as  defined  by  Suchanti  and  Fournier  for 
the  i'th  loop. 

the  error  integral  of  the  i'th  output  for  a  unit  step  change  in 

the  i'th  input  with  loop  i  closed  and  all  other  loops  open. 

the  error  integral  of  the  i'th  output  for  unit  step  changes  in 

all  inputs  with  all  loops  closed 

feedback  transfer  function  matrix. 

compensator  transfer  function  matrix. 

controller  transfer  function  matrix. 

parallel  transmittance  between  the  i'th  input  and  the  i'th  output 
when  the  i'th  loop  is  open  and  all  other  loops  are  closed 

steady  state  relative  gain  array, 
relative  dynamic  gain  array, 
average  dynamic  gain  array. 

matrix  used  in  definition  of  relative  transient  response 


functions. 
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Q={qijl 


R=hijt 

RL=frLij} 

R^=KiP 

S=diag{sj} 

u=[uii 

u(k) 

y^Cyil 

y(k) 

x(k) 

K 

*i 


CO 


compensated  plant  transfer  function  matrix. 

inverse  compensated  plant  transfer  function  matrix. 

minor  of  Q(s)«a,(s)  formed  by  deleting  all  rows  except  row  i 

and  all  columns  except  column  j. 

vector  of  reference  input  variables. 

closed  loop  transfer  function  matrix. 

closed  loop  load  transfer  function  matrix. 

transfer  function  matrix  for  system  with  all  but  one  feedback 
loop  closed. 

load  transfer  function  matrix  for  system  with  all  but  one 

feedback  loop  closed. 

feedback  TFM  with  s-=0  for  i=j,  and  Sj=  1  for  i^j. 
vector  of  compesated  plant  inputs, 
vector  of  inputs  to  discrete  state  space  model, 

vector  of  plant  output  variables. 

vector  of  outputs  from  discrete  state  space  model, 
vector  of  states  of  discrete  state  space  model, 

interaction  coefficient  as  defined  by  Rijnsdorp  (1965). 
extended  interaction  coefficient  for  the  i'th  input-output  pair 
of  a  m  x  m  multivariable  system, 
vector  of  load  or  disturbance  variables. 

the  relative  transient  response  function  betweein  the  j'th  input 

and  the  i'th  output. 

the  ratio  of  the  cofactor  of  element  (ij)  to  the  cofactor  of 

element  (ij)  of  the  return  difference  matrix, 
the  ratio  of  the  cofactor  of  element  (i,j)  to  the  cofactor  of 

element  (ij)  of  the  compensated  plant  transfer  function  matrix, 
the  continuous  frequency. 
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7.4  Nomenclature 

A=K1 

A 

C^Cijl 


D^diagJdjl 


for  chapter  3. 

real  matrix  with  positive  diagonal  elements  and  nonpositive 
offdiagonal  elements. 

coefficients  of  contraint  equations  of  general  linear 

programming  problem. 

companion  matrix  of  a  complex  matrix,  b;;=lq;;l  and  b-^=  —  I q j j I 
for  iAj. 

real  matrix  with  diagonal  elements  zero,  and  offdiagonal 
elements  nonnegative. 

element  of  diagonal  similarity  transformation  matrix,  which  will 
make  the  system  row  diagonal  dominant, 
diagonal  similarity  transformation  matrix. 


F 

g(x) 


G 


J 


N=hij! 


p. 

i 


Q=hijt 


xi 


W=diag^w;i 


return  difference  matrix. 

function  of  decision  variables  in  general  linear  programming 
problem. 

plant  transfer  function  matrix. 

transfer  function  between  the  i'th  input  and  the  i'th  output 
when  the  i'th  loop  is  open  and  all  other  loops  are  closed, 
function  to  be  maximized  in  transfer  of  dominance  or  general 
linear  programming  problem, 
polynomial  transfer  function  matrix. 

real  numbers,  which  are  indicators  of  column  matrix 
dominance. 

real  numbers,  which  are  indicatiors  of  row  matrix  dominance 
sum  of  off-diagonal  elements  of  row  i  of  a  matrix, 
complex  matrix. 

sum  of  off-diagonal  elements  of  column  i  of  a  matrix, 
vector  of  decision  variables  in  general  linear  programming 
problem. 

diagonal  matrix  with  positive  diagonal  elements. 

weighting  factor  for  i'th  row  or  i'th  column  in  transfer  of 
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dominance. 

V* 

function  obtained  by  mapping  of  the  Nyquist  O-contour  by 

q«- 

Q" 

normalized  minimal  Gershgorin  radius  as  calculated  by  Kantor 

and  Andres  method. 

CO 

continuous  frequency. 

7.5  Nomenclature 

for  chapter  4. 

matrix  of  derivatives  of  the  Nyquist  exact  loci  with  respet  to 

the  controller  gains. 

b=^b;} 

vector  of  the  distances  the  Nyquist  exact  loci  are  from  their 

desired  locations. 

G(s),G(z),or  G(w) 

plant  transfer  function  matrix. 

Gk(s) 

continuous  transfer  function  of  a  zero  order  hold. 

hi 

transfer  function  between  the  i'th  input  and  the  i'th  output  of 

a  multivariable  system  in  which  the  i'th  loop  is  open  and  all 

other  loops  are  closed. 

k*Di 

derivative  gain  in  the  i'th  loop. 

kIi 

integral  gain  in  the  i'th  loop. 

kp; 

proportional  gain  in  the  i'th  loop. 

•<1 

KA=diag^kjl 

compensator  transfer  function  matrix. 

controller  transfer  function  matrix. 

M 

Ikjhjl/1 1 +kjh-l  at  the  resonance  frequency. 

n ,  . 

hi 

number  of  clockwise  encirclements  of  the  critical  point  (-1,0) 

by  the  Nyquist  exact  loci  of  kjhj. 

Po 

number  of  poles  of  the  open  loop  transfer  function  matrix 

in  the  closed  right  half  plane. 

Q={q^ 

R6 

compensated  plant  transfer  function  matrix. 

band  aroung  the  Nyquist  loci  of  q  j;  within  which  the  loci  of 

the  transfer  function  h;  iies. 

V 

rise  time. 

rise  time. 
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tsi 

settling  time  to  within  1%  of  the  final  value 

T 

sampling  time. 

V 

discrete  frequency. 

x=ixi3 

vector  of  controller  gain  increments. 

*c 

relaxaton  factor  used  in  updating  the  controller  gains. 

i 

V,; 

overshoot. 

function  obtained  by  the  mapping  of  the  Nyquist  D~contour 

by  h,-. 

CO 

continuous  frequency. 

frequency  at  which  lk;h;l/l  1 +kjhjl  =  0.707,  the  bandwidth. 

<^>x 

frequency  at  which  lk,h;l/l  1  +kjh;l  =  0.5. 

cofactor  of  element  (i,j)  of  the  return  difference  matrix 

divided  by  the  cofactor  of  elememt  (i,i)  of  the  return 

difference  matrix. 

7.6  Nomenclature 

for  chapter  5. 

d 

characteristic  polynomial  of  transfer  function  matrix. 

F 

return  difference  matrix. 

FD0320 

constant  compensator/controller  designed  by  Kuon  (1975)  to 

give  a  gain  margin  of  5. 

FD0330 

constant  compensator/controller  designed  by  Kuon  to  give  gain 

margins  of  3.3. 

h  . 

i 

plant  transfer  function  matrix. 

transfer  function  between  the  i'th  input  and  the  i'th  output  of 

a  multivariable  system  in  which  the  i'th  loop  is  open  and  all 

other  loops  are  closed. 

Ki 

K*=diag{k;J 

constant  compensator  matrix. 

proportional  controller  gains. 

polynomial  part  of  a  transfer  function  matrix  with  a  common 

denominator  of  all  elements. 

n5. 

VJ 

real  numbers.  which  are  indicators  of  column  matrix 
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dominance. 

real  numbers,  which  are  indicators  of  row  matrix  dominance. 
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9.  Appendices 


9.1  Appendix  A.  Derivation  of  h  (s). 

In  this  appendix  expressions  for  the  transfer  function  hj(s)  between  the 

i'th  input  and  the  i'th  output  when  the  i'th  loop  is  open  and  all  other  loops  are 
closed  are  derived. 

Consider  the  multivariable  system  represented  by  the  block  diagram  in 
figure  9  1.  Q(s)Ka(s)  is  the  open  loop  transfer  function  matrix  and  S  is  a  dia¬ 
gonal  feedback  matrix,  defined  by  S  =  diag (s^] ,  k  =  1,...,m  with  Sj=0  and  s  =1, 

k£i,  i.e  the  i'th  feedback  loop  is  open.  The  closed  loop  transfer  function  for 

this  system  is 

H'(s)  =  O+Q'sIK^sISf1  Q(s)K-(s)  =  - - - fa'-te)}1’  Q(s)K*(s)  (9.1) 

det(l+Q(s)K^(s)S)  4 


where  det(  * )  denotes  the  determinant,  and  c~(s)  is  the  (i,j)'th  cofactor  of 
\+Q(s)K^(s)S.  The  sought  transfer  function  h((s)  is  the  (i,i)'th  element  of  H'(s), 
hence 


hj(s)  =  tys) 


1  m 

- T.  c..(s)q..(s) 

det(l+Q(s)KxS)  j=1  J 


(9.2  ) 


The  matrices  l+Q(s)Ka(s)S  and  l+Q(s)Ka(s)  only  differ  in  the  i'th  column.  This 
means  all  cofactors  of  the  i'th  column  of  l+Q(s)K^(s)S  are  the  same  as  the 
cofactors  of  the  i'th  column  of  l+CMsIK^is).  Denoting  the  cofactors  of  l+Q(s)K^(s) 
by  Cjjls)  (no  prime),  then 
m  m 

ZL  c;i(s)q:j(s)  =  ZI  c::(s)q  jj(s)  (9.3) 

j=1  J  J  j=1  J  J 


Furthermore,  Laplace  expansion  of  the  determinant  of  l+Q(s)K^(s)S  along  the  i'th 
column  gives 

det(l+Q(s)K*(s)S)  =  cH(s)  =  cu(s)  (94) 

since  only  the  i'th  element  of  the  i'th  column  of  l+Q(s)Ka(s)S  is  different  from 
zero.  This  element  is  one.  Hence  the  expression  for  hj(s)  becomes 
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Figure  9. 1 


Biockdiagram  of  multivariable  system  with  open  loop  transfer  function 
marix  Q(s)K  (s),  and  feedback  matrix  S  =  diag£svl,  s;=0,  sK=1,  k#i, 
i.e  all  loops  except  the  i'th  loop  are  closed. 
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h.(s) 


1 

C-.(s) 


m 


q^s) 


+ 


m 

X  4>ji  (s)q-j(s) 
j=  1  J  J 

tfi 


(9.5) 


The  above  expression  gives  h.(s)  in  terms  of  column  sums.  An  equivalent 
expression  can  be  established  in  terms  of  row  sums.  Laplace  expansion  of  the 
determinant  of  l+Qis)K^(s)  along  the  i'th  column  and  along  the  i'th  row  gives 

m  m 


det(l+Q(s)Ka(s))  = 


cii(s) 


j=i 


=  ci;(s) 


+ 


j=1 


Cj.(s)k;(s)qj;(s) 

M  J  *4 


(9.6) 


The  substitution  of  the  relationship 

m  m 

Cj.(s)kj(s)q..(s)  =  c.;(s)kj(s)qr( 

j=i  j1  1  y  j=  i  m  j  y 


(9.7) 


into  equation  9.5  gives  the  following  alternative  expression  for  h.(s): 


hj(s)  - 


1  m  m 

-  cjj(s)kj(s)q”(s)/k;(s)  =  q.js)  +  (J>-j <s)kj<s)qjj(s)/k|(s)  (9.8) 

Cii  s  J  ]*■< 


Equations  9.5  and  9.8  are  equivalent  to  the  inductively  derived  equations  (5.18) 
of  Kuon  (1975).  However,  the  normalization  of  the  return  difference  matrix  used 
by  Kuon  is  unnecessary,  and  in  equation  (5. 18e)  the  transpose  sign  should  be 
removed. 

Note  also,  the  above  hj(s)  is  different  from  the  hj(s)  introduced  by 
Rosenbrock  (1972).  Rosenbrock  used  hf(s)  as  the  transfer  function  between  the 
i'th  input  and  the  i'th  output  with  loops  1 i—  1  closed  and  loops  i . m  open. 
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9.2  Appendix  B.  Cofactor  equivalence. 

In  this  appendix  it  is  shown,  that  the  ratio  of  cofactors  of  the  return 

difference  matrix  F(s)  approach  the  corresponding  ratio  of  cofactors  of  the 

compensated  plant  transfer  function  matrix  Q(s)  in  the  limit  of  high  feedback 

gains. 

2x2  system:  Consider  the  matrices  Q(s)  and  F(s)  for  a  2  x  2  system 


Q(s)  = 

9ii 

9ia,~ 

F(s)  = 

1  +  Mu  k^q  1a_ 

9*a_ 

k-i^ai  1  +  k^q^ 

m 

with  the  following  ratios  of  cofactors  respectively 


9a,i  ^9 ** 

-a 

11 

"MM1  +  kaq^) 

(9.10) 

^1  = 

“9i*  /9n 

-MM'  +  k,q„) 

(9.1  1) 

Hence,  if  k  =  k^ 

=  kz,  and  k-*-»  then 

lim 

11 

lim 

<J>ai  =^*1 

(9.12) 

3x3  system:  Consider  the  ratios  of  cofactors  of  element  (1,2)  of  a  3  x  3 
system.  These  are  respectively 


9*3  9 31 

9  **9  33  “  P*S,C|3* 


(9.13) 


and 


k1k3<9^1933  q^3q3J  )  +  ^9^ 

(p^  _  -n  ”  — 

k*ka(qaaqs4  "  qaa^W  +  1  +  Ms*  +  Mas 

Hence,  if  k  =  k1  =  k^  =  k3,  and  k-*»«?  then 

iim  4>i*  =  'H'i* 
k  -*<*=> 


(9  14) 


(9.15) 


and  similarly  for  the  other  cofactor  ratios. 

m  x  m  system:  Let  Ka  =  kl,  and  k**«»  then  in  the  ratio  of  cofactors  of  I  + 

fvv.  ^  1 

klQ  only  those  elements  with  k  as  a  common  factor  will  not  vanish  in  the 
limit,  but  those  elements  are  exactly  the  ratio  of  cofactors  of  klQ  However, 


. 
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the  ratio  of  cofactors  of  klQ  and  Q  are  identical,  since  the  minors  of  kl 


cancel  out. 
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9.3  Appendix  C.  Proof  of  theorems  3.1  and  3.3. 

In  this  appendix  proofs  are  given  of  theorems  3.1  and  3.3  of  chapter  3. 

Proof  of  theorem  3.1:  Let  A  =  [a..}  be  the  companion  matrix  of  Q,  then 

•J 

according  to  the  definition  of  matrix  dominance  A  is  an  M-matrix,  and 
according  to  theorem  4.3  of  Fiedler  and  Ptak  (1962),  there  exist  a  diagonal 

matrix  D  =  diag{dj},  with  d;  >  0  for  all  i,  such  that  AD  is  row  diagonally 

dominant,  i.e  for  i  =  1,....,m 

m 

Iqi-.ldj  >  X  Iq-Jd.  (9.16) 

j=1  J  J 
J?i 

This  implies  that 

m 

Iqjjl  >  21  'qijKdj/dj)  (9.17) 

jTi 

The  lase  inequality  implies  QD  is  row  diagonally  dominant.  This  proves  the 

theorem. 

Proof  of  theorem  3.3:  Diagonal  dominance  of  D""1  QD  implies,  that  for  i  = 

1 . ,m 

m 

lq,*{l  >  !2L  IqjjKdj/d-)  (9.18) 

j=  1  J  J 

j^i 

for  column  diagonal  dominance,  and 

m 

Iq.-jl  >  51  lqr.Kd./d,)  (9. 19) 

j=1  J  J 

j*' 

for  row  diagonal  dominance.  These  are  equivalent  to  respectively 

m 

lq..|(1/d;)  -  21  Iq-.ld/ck)  >  0  (9.20) 

41  i  -  -|  V  J 

j* 


and 


. 
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m 

Iqjjldj  -  T.  Iqi5ld.  >  0  (9.21) 

j=  1  J  J 

j*i 

Now  let  xt  =  {d1, . d; . ,d^3  apd  y*  =  {l/d^, . 1/d; . ,  1/d,*}.  Then  inequalities 

9.20  and  9.21  can  be  rewritten  as: 

AT  y  >  0  (9.22) 

and 

A  x  >  0  (9.23) 

where  A  is  the  companion  matrix  of  Q.  Since  Q  is  row  (column)  diagonally 
dominant,  then  AT  (A)  is  an  M-matrix  (or  a  matrix  of  class  K),  and  a  solution  y 
(x)  to  inequalities  9.22  and  9.23  always  exist  (theorem  4.3  of  Fiedler  and  Ptak 
(1962)).  The  sought  matrix  D  is  then  the  matrix  with  the  elements  of  y  (x)  as 
diagonal  elements,  and  the  duality  theorem  is  proven. 
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9.4  Appendix  D.  Nyquist  contours  and  encirclements. 

In  this  appendix  the  modified  Nyquist  contours  in  the  s-plane,  z-plane 

and  w-plane  are  described  and  the  mapping  of  a  pure  integrator  is  described 
in  the  three  planes. 

The  Nyquist  contour  in  the  s-plane,  shown  in  figure  9.2  and  consisting 

of  the  imaginary  axis  and  a  large  semicircle  in  the  right  half  plane,  must  be 
modified  if  the  system  has  poles  on  the  imaginary  axis.  The  modified  Nyquist 

contour  in  the  s-plane  consist  of  a  large  semicircle  in  the  right  half  plane  and 

the  imaginary  axis  with  small  semicircles  around  the  poles  on  the  axis.  The 

small  semicircles  can  be  drawn  either  to  the  right  or  to  the  left  of  the 
imaginary  axis,  thus  giving  rise  to  the  two  different  modified  Nyquist  contours 

shown  in  figure  9.3.  The  contour  in  figure  9.3a  with  small  semicircles  to  the 

right  of  the  imaginary  axis  is  known  as  the  open  right  half  plane  Nyquist 

contour,  D^,  and  the  contour  in  figure  9.3b  with  small  semicircles  to  the  left 
of  the  imaginary  axis  is  known  as  the  closed  right  half  plane  Nyquist  contour, 
Dz.  The  modification  of  the  Nyquist  contour  is  necessary  for  the  contour  to 
remain  closed  and  the  mapping  of  the  contour  to  remain  conformal,  when  the 
mapping  function  has  poles  on  the  imaginary  axis.  Conformality  of  the  mapping 
is  necessary  to  allow  counting  of  encirclements. 

The  contours  in  the  z-plane,  D\  and  Dx,  are  shown  in  figure  9  4.  By 

the  transformation  z  =  es^*,  where  T  is  the  sample  time,  the  positive  imaginary 

axis  is  mapped  onto  the  upper  semicircle,  the  negative  imaginary  axis  is  mapped 
onto  the  lower  semicircle,  and  the  large  right  half  plane  semicircle  is  mapped 
onto  a  large  circle  in  the  z-plane.  The  contours  are  closed  by  a  cut  from 
(-1,0)  to  (-°°;0)  along  the  negative  real  axis.  Mathematically  the  large  semicircle  i 

the  s-plane  and  the  cut  and  the  large  circle  in  the  z-plane  represents  the  point 
at  infinity. 

The  Nyquist  contours  in  the  w-plane  ,  D*  and  Dx,  are  shown  in  figure 
9.5.  These  contours  represent  the  mapping  of  d\  and  by  the  bilinear  trans¬ 
formation  z  =  (1+w)/(1-w).  The  unit  circle  is  mapped  onto  the  imaginary  axis. 


. 
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A 


Figure  9.2  Nyquist  £>-contour  in  the  s-plane 


Figure  9.3  Modified  Nyquist  contour  in  the  s-plane.  a)  D*  -contour 
D-  -contour.  1 

'V 


■> 


and  b) 
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Figure  9.4  Modified  Nyquist  contours  in  the  z-plane.  a)  D'^  -contour,  and  b) 
D'  -contour. 


Figure  9.5  Modified  Nyquist  contours  in  the  w-plane.  a)  D*  -contour,  and  b) 
D£  -contour. 
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the  (-1,0)  point  is  mapped  onto  the  point  at  infinity,  the  cut  along  the  negative 
real  axis  is  mapped  onto  a  cut  along  the  positive  real  axis  from  («*?0)  to  (1,0). 
The  contours  are  closed  by  two  quartercircles  representing  the  point  at  infinity. 

The  use  of  the  different  modified  Nyquist  contours  is  illustrated  by 
considering  a  pure  integrator.  The  mappings  of  DA  and  D%  by  the  continuous 
transfer  function  G(s)  =  1/s  are  shown  in  figure  9.6.  As  expected  the  mapping 
of  gives  zero  encirclements  and  the  mapping  of  Dx  gives  one  encirclement 
of  tne  origin  and  the  point  (-1,0).  The  discrete  integrator,  or  the  z-transform 
of  GjJs)G(s)  =  (1-esT)/s2',  is  G(z)  =  T/(z-1).  The  mappings  of  D\  and  D x  are 
shown  in  figure  9.7.  Again  the  mapping  of  D\  gives  zero  encirclements  and  the 
mapping  of  Dx  gives  one  encirclement  of  the  origin  and  the  point  (-1,0). 
However,  this  result  hinges  on  the  careful  mapping  of  the  cut  along  the 
negative  real  axis.  Offen  G(z)  is  subjected  to  a  bilinear  transformation  z  = 
( 1  +  w)/ ( 1  —  w)  before  mapping  of  the  Nyquist  contour.  With  this  transformation 
G(z)  becomes  G(w)  =  -T/2  +  T/2w.  The  mappings  of  D\  and  D ^  by  G(w)  are 
shown  in  figure  9.8.  These  are  identical  to  those  in  figure  9.6,  due  to  the 
mapping  of  the  cut  along  the  positive  real  axis.  This  shows,  it  is  crucial  for 
the  correct  encirclement  result  to  be  obtained,  that  the  contours  d\  or  D*  be 
used,  and  not  D1  or  Dx,  when  bilinear  transformed  systems  are  considered. 
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Figure  9.6  Mappings  of  D  A  and  Dx  by  G(s)  =  1/s. 


-> 


-> 


Figure  9  7  Mappings  of  D'  and  D'  by  G(z)  =  T /(z—  1 ). 


'M  * 


and  d\  by  G(w)  =  -T/2  +  T/2w. 


Figure  9.8  Mappings  of  D\ 
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9.5  Appendix  E.  Proof  of  theorem  4.1. 

In  this  appendix  a  proof  is  given  of  theorem  4. 1  of  chapter  4. 

Proof  of  theorem  4.1:  Let  cjj(s)  be  the  cofactor  of  element  (i,j)  of  the  return 

>J 

difference  matrix  F(s).  From  equation  9.5  of  appendix  A 

IF(s)l 

1  +k.(s)h.(s)  =  -  (9.24) 

c-(s) 

n 

hence 

IF(s)l  (9.25) 

Since  F(s)  is  matrix  dominant,  by  theorem  3.1  the  mapping  of  the  Nyquist 

£)-contour  by  the  determinant  IF(s)l  gives  the  same  number  of  encirclements  of 

the  origin  as  the  sum  of  encirclements  in  the  mapping  of  D  by  the  diagonal 
elements  fj,*(s)  of  F(s).  Let  the  mapping  of  D  by  f  (s)  encircle  the  origin  n j j 

times  clockwise,  then  the  mapping  of  D  by  IF(s)l  will  give  rise  to  nT1 
clockwise  encirclements  of  the  origin,  with 

m 

nT<1  =  (m-1)  nj;  (9.26) 

i=  1 

Since  the  minors  of  the  diagonal  elements  of  the  matrix  F(s)  and  the  diagonally 
dominant  matrix  D(s)-i  F(s)D(s)  (with  D(s)  =  diag{d-((s)},  dj(s)  >  0)  are  the  same, 
the  mapping  of  D  by  c--(s)  gives  the  following  number  of  clockwise 

encirclements  of  the  origin 
m 

m.  =  T.  n;i  (9.27) 

j  j=i  " 

i^j 

m. 

So  the  mapping  of  D  by  1?  c-(s)  will  give  the  following  number  of 

i  =  1 


m 


TT  (1+k,(s)h;(s))  = 
i=  1 


IF(s)l 


-1 


m 

7T 

i=l 


c-  s 

it 


encirclements: 


. 
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m  m 


m 


n 


rz 


(m-1)  H  n  - 


(9.28) 


Therefore  n.^.  =  n-p^,  and  by  the  principle  of  the  argument  the  net  number  of 
encirclements  by  the  quantity  in  the  square  bracket  in  equation  9.25  is  zero. 
Hence  the  number  of  encirclements  of  the  critical  point  (-1,0)  by  the  mapping 

of  D  by  kj(s)hj(s)  ,  i  =  1, m  is  equal  to  the  number  of  encirclements  of  the 

origin  by  the  mapping  of  D  by  IF(s)l,  provided  F(s)  is  matrix  dominant.  This 
proves  the  theorem. 


